Math 4401 Numerical Methods

Notebook Author/Course Instructor: Barry McQuarrie
Date: Sept 2016
Lecture topic: Newton's Method for Systems

Nonlinear Systems

An Example of a Nonlinear System in (R^2)

Here is a nonlinear system of equations

f[x_,y_]=x^2y-y^3-1/2
g[x_,y_]=Cos[x]y x-y^2x-1

Here are plots of the two equations in blue and red. Where the two curves intersect, we have a solution to the system of equations.

plot1=ContourPlot[f[x,y]==0,{x,-20,20},{y,-10,10},ContourStyle->{Blue,Thick},PlotPoints->300];
plot2=ContourPlot[g[x,y]==0,{x,-20,20},{y,-10,10},ContourStyle->{Red,Thick},PlotPoints->300];
Show[plot1,plot2,PlotRange->{{-10,10},{-10,10}},ImageSize->{500,Automatic},BaseStyle->20,AxesLabel->{"x", "y"}]

The intersection of the lines will be the solutions to the system. Very difficult to find! Looks from our picture that we have two real valued solutions. This would be difficult to do in higher dimensions.

Mathematica has difficulty with these equations:

NSolve[{f[x,y]==0,g[x,y]==0},{x,y}]

By the time we are done, we will have a solution.

Norm, Guassian elimination module, and an inverse matrix module

We will need to take the inverse of a matrix, so that means a Guassian Elimination module and inverse matrix module are needed. We will also need a norm, to be able to tell how our vector x is converging.

norm2[g_]:=((((∑_(i=1)^(Length[g]))((g[[i]])^2)))^(1/2));

GaussElim3[A0_]:=Module[ {A=A0,p,k,i,n=Length[A0]},  
For[ p=1, p≤n, p=p+1, 
For[ k=p+1, k≤n, k=k+1, 
If[ Abs[A[[k,p]] ]>Abs[A[[p,p]] ], A[[{p,k}]]  = A[[{k,p}]] ;   ]; ];
A[[p]] = ((A[[p]])/(A[[p,p]] )); 
For[ i=1, i≤n, i=i+1, If[ i≠p, A[[i]] = A[[i]] - A[[i,p]] A[[p]]  ]; ]; 
 ]; 
Return[A];  ] ;

MatInv[A0_]:=Module[{A1,A2,A3,n=Length[A0]},
A1=Join[A0,IdentityMatrix[n],2];
A2=GaussElim3[A1];
A3=Take[A2,{1,n},{n+1,2n}];
Return[A3]
];

The Jacobian (the gradient vector field)

The Jacobian matrix is the matrix of partial derivatives of a function.

It can be defined in Mathematica using the Outer command.

A0 is the list of equations, and vars is the list of variables.

jacobian[A0_,vars_]:=Outer[D,A0,vars]

Here is the general case:

jacobian[{F[x1,x2,x3],G[x1,x2,x3],H[x1,x2,x3]},{x1,x2,x3}];
MatrixForm[%]

Here is the Jacobian for our equations:

F={f[x,y],g[x,y]}
var={x,y}

jac[x_,y_]=jacobian[F,var];
MatrixForm[%]

At a point in space, the Jacobian tells you what the functions would like to do.

The rows of the Jacobian are the components of the gradient vector field for the function.

vecplot1=VectorPlot[{jac[x,y][[1,1]],jac[x,y][[1,2]]},{x,-20,20},{y,-10,10},VectorPoints->40,VectorScale->{0.02,Large,None}];
Show[plot1,vecplot1,PlotRange->{{-5,10},{-5,10}},Frame->True,ImageSize->{500,Automatic}]

vecplot2=VectorPlot[{jac[x,y][[2,1]],jac[x,y][[2,2]]},{x,-20,20},{y,-10,10},VectorPoints->40,VectorScale->{0.02,Large,None}];
Show[plot2,vecplot2,PlotRange->{{-5,10},{-5,10}},Frame->True,ImageSize->{500,Automatic}]

We can zoom in to see what is happening around the "circular" region.

vecplot2b=VectorPlot[{jac[x,y][[2,1]],jac[x,y][[2,2]]},{x,5,7},{y,0,1},VectorPoints->20,VectorScale->{0.02,Large,None}];
plot2b=ContourPlot[g[x,y]==0,{x,5,7},{y,0,1},ContourStyle->{Red,Thick},PlotPoints->300];
Show[plot2b,vecplot2b,PlotRange->{{5,7},{0,1}},Frame->True,ImageSize->{700,Automatic},AspectRatio->1/2]

In all these plots, we notice that the vector <(Subscript[f, x])(x,y),(Subscript[f, y])(x,y)> intersects the curve F(x,y)=0 at right angles.

This vector, the gradient vector, gives the direction of fastest increase in F.

The curve F(x,y)=0 is a level curve of F, and the gradient vector intersects the level curves orthogonally.

The Vector Field - (([J(x) ])^(-1))F(x)

Here are the functions we want to solve, and the Jacobian:

F
jac[x,y]

Here is the vector field that enters Newton's method for systems:

-MatInv[jac[x,y]].F//Simplify//MatrixForm

Here is the vector field plot of the above quantity:

vecplote=VectorPlot[-MatInv[jac[x,y]].F,{x,-10,5},{y,-5,5},VectorPoints->30,VectorScale->{0.02,Large,None}];
Show[plot1,plot2,vecplote,PlotRange->{{-5,5},{-5,5}},Axes->False,AxesOrigin->{0,0},Axes->Automatic,AspectRatio->Automatic,ImageSize->{700,Automatic}]

Zooming in will help clarify some regions:

vecplotd=VectorPlot[-MatInv[jac[x,y]].F,{x,-1,0},{y,-1.5,-0.5},VectorPoints->30,VectorScale->{0.02,Large,None}];
Show[plot1,plot2,vecplotd,PlotRange->{{-1,0},{-1.5,-0.5}},Axes->False,AxesOrigin->{0,0},Axes->Automatic,AspectRatio->Automatic,ImageSize->{700,Automatic}]

Notice that the vector field - (([J(x) ])^(-1))F(x) converges at the point of intersection we seek! This is truly excellent.

vecplotd=VectorPlot[-MatInv[jac[x,y]].F,{x,-2,-0.1},{y,0.1,2},VectorPoints->30,VectorScale->{0.02,Large,None}];
Show[plot1,plot2,vecplotd,PlotRange->{{-2,0},{0,2}},Axes->False,AxesOrigin->{0,0},Axes->Automatic,AspectRatio->Automatic,ImageSize->{700,Automatic}]

We can have problems if the vector field requires the inverse of a matrix which has no inverse.

This is the equivalent of division by zero when the derivative equals zero in the one dimensional case.

vecplotf=VectorPlot[-MatInv[jac[x,y]].F,{x,-10,10},{y,-10,10},VectorPoints->30,VectorScale->{0.02,Large,None}];
Show[plot1,plot2,vecplotf,PlotRange->{{-10,10},{-10,10}},Axes->False,AxesOrigin->{0,0},Axes->Automatic,AspectRatio->Automatic,ImageSize->{700,Automatic}]

If you begin close to a root, the vector field will push you towards the solution.

If you do not begin close to a root, the vector field can bounce you around quite dramatically!

Newton's Method for Nonlinear Systems

Newton's Method: Graphical Output

Define the set of functions we wish to find the zero of:

f1[x_,y_]=F

And here is Newton's Method for systems. The only change from before is using an inverse matrix module, and using matrix multiplication.

Clear[p0,p1,tol, max, i,p,q,y,yn]

p0={-10.0,-3.0};
tol=10.0^(-10);
max=20;

ptsize=0.08
data={Graphics[{Darker[Red],Disk[{p0[[1]],p0[[2]]},ptsize]}]};
data2={p0};
k=1;
Print["Newton's Method for Systems"]
Print["Input:  ","p0= ",p0," tol= ",tol]
While[k≤max ,

yn= MatInv[jac[p0[[1]],p0[[2]]]].f1[p0[[1]],p0[[2]]];
p = p0 -yn;

data=Append[data,Graphics[{Darker[Red],Disk[{p[[1]],p[[2]]},ptsize]}]];
data2=Append[data2,p];
If[norm2[yn]≤tol,{Print["Result: ","p= ",SetPrecision[p,10]," Error ",norm2[p-p0], " Steps= ",k],Break[]}];
p0=p;
k=k+1]
If[k>max,Print["Method Failed to Converge"]]

Show[plot1,plot2,vecplotf,Graphics[{Darker[Green],Thick,Line[data2]}],data,PlotRange->{{-10,0},{-10,0}},Axes->True,AxesOrigin->{0,0},Axes->Automatic,AspectRatio->Automatic,ImageSize->{500,Automatic}]

The red points are the output from Newton's Method. The green line shows how Newton's method moves from point to point, based on a "kick" from the vector field.

Newton's Method: Graphical Output

Define the set of functions we wish to find the zero of:

f1[x_,y_]=F

And here is Newton's Method for systems:

Clear[p0,p1,tol, max, i,p,q,y,yn]

p0={-9.0,-2.1};
tol=10.0^(-10);
max=100;

ptsize=0.08
data={Graphics[{Darker[Red],Disk[{p0[[1]],p0[[2]]},ptsize]}]};
data2={p0};
k=1;
Print["Newton's Method for Systems"]
Print["Input:  ","p0= ",p0," tol= ",tol]
While[k≤max ,

yn= MatInv[jac[p0[[1]],p0[[2]]]].f1[p0[[1]],p0[[2]]];
p = p0 -yn;

data=Append[data,Graphics[{Darker[Red],Disk[{p[[1]],p[[2]]},ptsize]}]];
data2=Append[data2,p];
If[norm2[yn]≤tol,{Print["Result: ","p= ",SetPrecision[p,10]," Error ",norm2[p-p0], " Steps= ",k],Break[]}];
p0=p;
k=k+1]
If[k>max,Print["Method Failed to Converge"]]

Show[plot1,plot2,vecplotf,Graphics[{Darker[Green],Thick,Line[data2]}],data,PlotRange->{{-10,1},{-5,6}},Axes->True,AxesOrigin->{0,0},Axes->Automatic,AspectRatio->Automatic,ImageSize->{500,Automatic}]

The red points are the output from Newton's Method. The green line shows how Newton's method moves from point to point, based on a "kick" from the vector field. Notice that in this case, we jumped around quite a bit before finding a solution to the system of equations.

Newton's Method: Table Output

Clear[p0,p1,tol, max, i,p,q]

f1[x_,y_]={-(1/2)+(x^2) y-(y^3),-1-x (y^2)+x y Cos[x]}

p0={-1.0,-1.0};
tol=10.0^(-10);
max=20;

data={};
k=1;
Print["Newton's Method for Systems"]
Print["Input:  ","p0= ",p0," tol= ",tol]
While[k≤max ,

y= MatInv[jac[p0[[1]],p0[[2]]]].f1[p0[[1]],p0[[2]]];
p = p0 -y;

data=Append[data,{k,p,norm2[y]}];
If[norm2[y]≤tol,{Print["Result: ","p= ",SetPrecision[p,10]," Error ",norm2[p-p0], " Steps= ",k],Break[]}];
p0=p;
k=k+1]
If[k>max,Print["Method Failed to Converge"]]

TableForm[SetPrecision[data,10],TableHeadings->{None,{"n","p","Norm p"}}]

Newton's Method for Nonlinear Systems Higher Dimension and Finding Good Starting Point

Quasi-Random: Corput Sequence in (ℝ^1)

Recall the base expansion for a number:

928.64=9×(10^2)+2×(10^1)+8×(10^0)+6×(10^(-1))+4×(10^(-2))

Here is a Mathematica commands you may not have seen before, which help us pull out the coefficients:

IntegerDigits[928,10]

This works in any base:

IntegerDigits[6,2]

The Corput sequence is defined as follows:

The decimal number n is expanded in base b.

The number in base b is reflected.

Convert back to decimal number.

where we do this process for n=1,2,3,4,... to generate the Corput sequence. Here is the Corput number calculated when n=6 and b=2:

Expand in base 2: 6=1×(2^2)+1×(2^1)+0×(2^0)-> 110 in base 2
Reflect in base 2: 011-> 0×(2^(-1))+1×(2^(-2))+1×(2^(-3))
Convert to decimal: (1/2)+(1/8)=(3/8)

This module automates the process outlined by hand:

Corput[n_,b_]:=Module[{expandbaseb,m},
expandbaseb=IntegerDigits[n,b];
m=Length[expandbaseb];
expandbaseb=Reverse[expandbaseb];
Return[(∑_(j=0)^(m-1))expandbaseb[[j+1]]*(b^(-j-1))];
]

Corput[6,2]

And here is the Corput sequence when the base is b=2.

corput=Table[Corput[i,2],{i,1,20}]

Here is a plot of how the Corput sequence with base 2 fills the region x ϵ (0,1). Notice that any “holes” in the region are filled by the sequence, which is what makes it a quasi-random sequence (it avoids clumping).

Clear[list]
n=80;

list={};
Do[
list=Append[list,Graphics[Disk[{Corput[i,2],0},{0.01,0.01}]]]
,{i,1,n}];

Manipulate[
Show[Take[list,i],PlotRange->{{0,1},{-(1/8),(1/8)}},Axes->{True,False},ImageSize->{500,Automatic},AspectRatio->(1/4)]
,{i,1,n,1}]

Here is a plot of how the Corput sequence with base 3 fills the region x ϵ (0,1).

corput=Table[Corput[i,3],{i,1,20}]

Clear[list]
n=80;

list={};
Do[
list=Append[list,Graphics[Disk[{Corput[i,3],0},{0.01,0.01}]]]
,{i,1,n}];

Manipulate[
Show[Take[list,i],PlotRange->{{0,1},{-(1/8),(1/8)}},Axes->{True,False},ImageSize->{500,Automatic},AspectRatio->(1/4)]
,{i,1,n,1}]

Quasi-Random: Halton Sequence in (ℝ^2)

To extend to higher dimensions, you take as base the next prime number, so in (ℝ^2) you use the x-coordinate with base 2, the y-coordinate with base 3.

Halton2D[n_]:={Corput[n,2],Corput[n,3]}

halton=Table[Halton2D[k],{k,1,10}]

Here are some examples of how the Halton sequence fills space in (ℝ^2).

halton=Table[Halton2D[k],{k,1,2000}];

hplot=ListPlot[halton,Frame->True,AspectRatio->1,PlotStyle->Red,ImageSize->{500,Automatic}]

Here is an animation of how the Halton sequence fills the x y-plane.

Clear[list]
n=1200;

list={};
Do[
list=Append[list,Graphics[{Red,Disk[Halton2D[i],{0.01,0.01}]}]]
,{i,1,n}];

Manipulate[
Show[Take[list,i],PlotRange->{{0,1},{0,1}},Frame->True,AspectRatio->1,ImageSize->{500,Automatic}]
,{i,1,n,1}]

Newton's Method:  Higher Dimension

Real Roots

In higher dimensions, in can be difficult to get good starting values. It can also occur that the Jacobain matrix does not have an inverse at some points.

Here is the system we want to solve

Clear[f1,x,y,z,jac,yn]
f1[x_,y_,z_]={10x-2y^2+y-2z-5,x+8y^2+4z^2-9,8y*z^2+4x}

jac[x_,y_,z_]=jacobian[f1[x,y,z],{x,y,z}]

We need to do some preprocessing to find a good starting point. Sample a small region of space where you think a solution lies.To find good starting values in general, you will need to do something more powerful, like use Steepest Descent. Here, we can use the quasi-random Corput sequence to sample space, and find the smallest norm at a variety of points. We will examine this Corput sequence when we look at pseudo-random numbers later in the semester, so seeing it here is just a little foreshadowing.

Corput[n_,b_]:=Module[{expandbaseb,m},
expandbaseb=IntegerDigits[n,b];
m=Length[expandbaseb];
expandbaseb=Reverse[expandbaseb];
Return[(∑_(j=0)^(m-1))expandbaseb[[j+1]]*(b^(-j-1))];
]

x1=-1.0;
x2=1.0;
y1=-1.0;
y2=1.0;
z1=-1.0;
z2=1.0;
minpts = {0,0,0};
min = 10.0;
Do[
{xp =x1+(x2-x1)*Corput[i,2]//N;
yp =y1+(y2-y1)*Corput[i,3]//N;
zp =z1+(z2-z1)*Corput[i,5]//N;
fcn1 =norm2[f1[xp,yp,zp]];
If[fcn1<min,{min=fcn1,minpts={xp,yp,zp}}]
},{i,1,1000}]

minpts

We will use this to initialize and find a solution to the equations.

Clear[p0,p1,tol, max, i,p,q]

p0=minpts;
tol=10.0^(-14);
max=20;

data={};
k=1;
Print["Newton's Method for Systems"]
Print["Input:  ","p0= ",p0," tol= ",tol]
While[k≤max ,

yn= MatInv[jac[p0[[1]],p0[[2]],p0[[3]]]].f1[p0[[1]],p0[[2]],p0[[3]]];
p = p0 -yn;

data=Append[data,{k,p,norm2[yn]}];
If[norm2[yn]≤tol,{Print["Result: ","p= ",SetPrecision[p,10]," Error ",norm2[p-p0], " Steps= ",k],Break[]}];
p0=p;
k=k+1]
If[k>max,Print["Method Failed to Converge"]]

TableForm[SetPrecision[data,10],TableHeadings->{None,{"n","p","Norm p"}}]

Let's check this against Mathematica's result:

NSolve[f1[x,y,z],{x,y,z}]

Complex Roots

We can also use complex numbers to initialize this procedure, and find complex roots. We will want to make sure our norm produces real valued numbers, so let's add an absolute value to it:

norm2[g_]:=((((∑_(i=1)^(Length[g]))((Abs[g[[i]]])^2)))^(1/2));

x1=0.0;
x2=1.0;
y1=0.0;
y2=1.0;
z1=0.0;
z2=1.0;
minpts = {0,0,0};
min = 10.0;
Do[
{xp =x1+(x2-x1)*Corput[i,2]+I*(x1+(x2-x1)*Corput[i,7])//N;
yp =y1+(y2-y1)*Corput[i,3]+I*(y1+(y2-y1)*Corput[i,11])//N;
zp =z1+(z2-z1)*Corput[i,5]+I*(z1+(z2-z1)*Corput[i,13])//N;
fcn1 =norm2[f1[xp,yp,zp]];
If[fcn1<min,{min=fcn1,minpts={xp,yp,zp}}]
},{i,1,1000}]

minpts

Clear[p0,p1,tol, max, i,p,q]

p0=minpts;
tol=10.0^(-14);
max=20;

data={};
k=1;
Print["Newton's Method for Systems"]
Print["Input:  ","p0= ",p0," tol= ",tol]
While[k≤max ,

yn= MatInv[jac[p0[[1]],p0[[2]],p0[[3]]]].f1[p0[[1]],p0[[2]],p0[[3]]];
p = p0 -yn;

data=Append[data,{k,p,norm2[yn]}];
If[norm2[yn]≤tol,{Print["Result: ","p= ",SetPrecision[p,10]," Error ",norm2[p-p0], " Steps= ",k],Break[]}];
p0=p;
k=k+1]
If[k>max,Print["Method Failed to Converge"]]

TableForm[SetPrecision[data,10],TableHeadings->{None,{"n","p","Norm p"}}]

Let's check this against Mathematica's result:

NSolve[f1[x,y,z],{x,y,z}]

Well done indeed. If you want to find other roots to the system, you would need to find different initializing values.

Fun Graphics: Solutions of Nonlinear (R^2) Systems with Real Roots 
(takes a while to run!)

Initializations

First, ensure some modules are loaded (a norm, a Gaussian Elimination, a matrix inverse, and the jacobian):

norm2[g_]:=((((∑_(i=1)^(Length[g]))((Abs[g[[i]]])^2)))^(1/2));

GaussElim3[A0_]:=Module[ {A=A0,p,k,i,n=Length[A0]},  
For[ p=1, p≤n, p=p+1, 
For[ k=p+1, k≤n, k=k+1, 
If[ Abs[A[[k,p]] ]>Abs[A[[p,p]] ], A[[{p,k}]]  = A[[{k,p}]] ;   ]; ];
A[[p]] = ((A[[p]])/(A[[p,p]] )); 
For[ i=1, i≤n, i=i+1, If[ i≠p, A[[i]] = A[[i]] - A[[i,p]] A[[p]]  ]; ]; 
 ]; 
Return[A];  ] ;

MatInv[A0_]:=Module[{A1,A2,A3,n=Length[A0]},A1=Join[A0,IdentityMatrix[n],2];
A2=GaussElim3[A1];A3=Take[A2,{1,n},{n+1,2n}];Return[A3]];

Corput[n_,b_]:=Module[{expandbaseb,m},
expandbaseb=IntegerDigits[n,b];
m=Length[expandbaseb];
expandbaseb=Reverse[expandbaseb];
Return[(∑_(j=0)^(m-1))expandbaseb[[j+1]]*(b^(-j-1))];
]

jacobian[A0_,vars_]:=Outer[D,A0,vars]

Zoomed Out

Here is where we define the equations to be solved (the same equations as before)

f1[x_,y_]={x^2y-y^3-1/2,Cos[x]y x-y^2x-1};
jac[x_,y_]=jacobian[f1[x,y],{x,y}];

Let's finish with a plot of the space that is used to initialize a system in (R^2), and which root Newton's method for systems finds.

Clear[acc,tol,max,nr,p0,p,set,roots]
acc=5;
tol=(10.^(-acc));
max=100;
nr=0;
x1=-5.;
x2=6.;
y1=-5.;
y2=5.;
numrun=5000;
roots={};
Do[set[i]={},{i,1,numrun}]; (*set[i] is the set of all initial points that found root[i]*)
Do[
xp=N[x1+(x2-x1) Corput[i,2]];
yp=N[y1+(y2-y1) Corput[i,3]];
p0={xp,yp};
data={};
k=1;
While[k≤max,
yn=MatInv[jac[p0[[1]],p0[[2]]]].f1[p0[[1]],p0[[2]]];
p=p0-yn;
If[norm2[yn]≤tol,
newroot=True;
Do[If[norm2[p0-roots[[i]]]<(10^(-acc+1)),newroot=False];,{i,1,Length[roots]}];
If[newroot,roots=Append[roots,p0]];Do[If[norm2[p0-roots[[i]]]<(10^(-acc+1)),set[i]=Append[set[i],{xp,yp}]],{i,1,Length[roots]}];
Goto[donethis];];
p0=p;
k=k+1];
nr=nr+1;
Label[donethis];,
{i,1,numrun}]
Print["Did not find a root for ",nr," initial conditions"] 

setlist={};
Do[setlist=Append[setlist,set[i]],{i,1,Length[roots]}] ;
huelist={};
Do[huelist=Append[huelist,Hue[(i/(Length[roots]))]],{i,1,Length[roots]}];phasepts=ListPlot[setlist,AspectRatio->1,Axes->True,Frame->True,PlotRange->{{x1,x2},{y1,y2}},PlotStyle->huelist,Background->GrayLevel[1],ImageSize->{600,Automatic}]

Here are plots of the two equations we are trying to solve in red and blue. Where the two curves intersect, we have a solution to the system of equations.

plot1=ContourPlot[f1[x,y][[1]]==0,{x,x1,x2},{y,y1,y2},ContourStyle->{Blue,Thick},PlotPoints->300];
plot2=ContourPlot[f1[x,y][[2]]==0,{x,x1,x2},{y,y1,y2},ContourStyle->{Red,Thick},PlotPoints->300];
vecplotf=VectorPlot[-MatInv[jac[x,y]].f1[x,y],{x,x1,x2},{y,y1,y2},VectorPoints->30,VectorScale->{0.02,Small,None}];
Show[phasepts,plot1,plot2,vecplotf,PlotRange->{{x1,x2},{y1,y2}},Axes->False,AxesOrigin->{0,0},Axes->Automatic,AspectRatio->1,ImageSize->{600,Automatic}]

Some parts of the sample space didn't find a root, and others did.

More Points (10 minutes in office)

Here we just sample more points:

Clear[acc,tol,max,nr,p0,p,set,roots]
acc=5;
tol=(10.^(-acc));
max=100;
nr=0;
x1=-5.;
x2=6.;
y1=-5.;
y2=5.;
numrun=40000;
roots={};
Do[set[i]={},{i,1,numrun}]; (*set[i] is the set of all initial points that found root[i]*)
Do[
xp=N[x1+(x2-x1) Corput[i,2]];
yp=N[y1+(y2-y1) Corput[i,3]];
p0={xp,yp};
data={};
k=1;
While[k≤max,
yn=MatInv[jac[p0[[1]],p0[[2]]]].f1[p0[[1]],p0[[2]]];
p=p0-yn;
If[norm2[yn]≤tol,
newroot=True;
Do[If[norm2[p0-roots[[i]]]<(10^(-acc+1)),newroot=False];,{i,1,Length[roots]}];
If[newroot,roots=Append[roots,p0]];Do[If[norm2[p0-roots[[i]]]<(10^(-acc+1)),set[i]=Append[set[i],{xp,yp}]],{i,1,Length[roots]}];
Goto[donethis];];
p0=p;
k=k+1];
nr=nr+1;
Label[donethis];,
{i,1,numrun}]
Print["Did not find a root for ",nr," initial conditions"] 

setlist={};
Do[setlist=Append[setlist,set[i]],{i,1,Length[roots]}] ;
huelist={};
Do[huelist=Append[huelist,Hue[(i/(Length[roots]))]],{i,1,Length[roots]}];phasepts=ListPlot[setlist,AspectRatio->1,Axes->True,Frame->True,PlotRange->{{x1,x2},{y1,y2}},PlotStyle->huelist,Background->GrayLevel[1],ImageSize->{600,Automatic}]

Here are plots of the two equations we are trying to solve in blue and red. Where the two curves intersect, we have a solution to the system of equations.

plot1=ContourPlot[f1[x,y][[1]]==0,{x,x1,x2},{y,y1,y2},ContourStyle->{Blue,Thick},PlotPoints->300];
plot2=ContourPlot[f1[x,y][[2]]==0,{x,x1,x2},{y,y1,y2},ContourStyle->{Red,Thick},PlotPoints->300];
vecplotf=VectorPlot[-MatInv[jac[x,y]].f1[x,y],{x,x1,x2},{y,y1,y2},VectorPoints->30,VectorScale->{0.02,Small,None}];
Show[phasepts,plot1,plot2,vecplotf,PlotRange->{{x1,x2},{y1,y2}},Axes->False,AxesOrigin->{0,0},Axes->Automatic,AspectRatio->1,ImageSize->{600,Automatic}]

Another System

Let's look at this system:

f1[x_,y_]={x^2/6+y^2/8-1,Cos[y]x+Cos[x](y^2)};
jac[x_,y_]=jacobian[f1[x,y],{x,y}];

x1=-5;x2=5;y1=-5.;y2=5.;

plot1=ContourPlot[f1[x,y][[1]]==0,{x,x1,x2},{y,y1,y2},ContourStyle->{Blue,Thickness[0.01]},PlotPoints->300];
plot2=ContourPlot[f1[x,y][[2]]==0,{x,x1,x2},{y,y1,y2},ContourStyle->{Darker[Red],Thickness[0.01]},PlotPoints->300];
vecplotf=VectorPlot[-MatInv[jac[x,y]].f1[x,y],{x,x1,x2},{y,y1,y2},VectorPoints->20,VectorScale->{0.03,Medium,None},VectorStyle->Thick];

Show[plot1,plot2,vecplotf]

NSolve[{f1[x,y][[1]]==0,f1[x,y][[2]]==0},{x,y}]

To speed things up a bit, let’s work out the jaobian bit in general so we aren’t inverting a matrix a ton of times:

jacbit[x_,y_]=MatInv[jac[x,y]].f1[x,y]

Here we zoom in on an interesting area:

Clear[acc,tol,max,nr,p0,p,set,roots]
acc=7;
tol=(10.^(-acc));
max=200;
nr=0;
x1=-5;
x2=5;
y1=-5.;
y2=5.;
numrun=40000;
roots={};
Do[set[i]={},{i,1,numrun}]; (*set[i] is the set of all initial points that found root[i]*)
Do[
xp=N[x1+(x2-x1) Corput[i,2]];
yp=N[y1+(y2-y1) Corput[i,3]];
p0={xp,yp};
data={};
k=1;
While[k≤max,
yn=jacbit[p0[[1]],p0[[2]]];
p=p0-yn;
If[norm2[yn]≤tol,
newroot=True;
Do[If[norm2[p0-roots[[i]]]<(10^(-acc+3)),newroot=False];,{i,1,Length[roots]}];
If[newroot,roots=Append[roots,p0]];Do[If[norm2[p0-roots[[i]]]<(10^(-acc+3)),set[i]=Append[set[i],{xp,yp}]],{i,1,Length[roots]}];
Goto[donethis];];
p0=p;
k=k+1];
nr=nr+1;
Label[donethis];,
{i,1,numrun}]
Print["Did not find a root for ",nr," initial conditions"] 

setlist={};
Do[setlist=Append[setlist,set[i]],{i,1,Length[roots]}] ;
huelist={};
Do[huelist=Append[huelist,Hue[(i/(Length[roots]))]],{i,1,Length[roots]}];phasepts=ListPlot[setlist,AspectRatio->1,Axes->True,Frame->True,PlotRange->{{x1,x2},{y1,y2}},PlotStyle->huelist,Background->GrayLevel[1],ImageSize->{600,Automatic}]

Here are plots of the two equations we are trying to solve in red and blue. Where the two curves intersect, we have a solution to the system of equations.

Show[phasepts,plot1,plot2,vecplotf,PlotRange->{{x1,x2},{y1,y2}},Axes->False,AxesOrigin->{0,0},Axes->Automatic,AspectRatio->1,ImageSize->{600,Automatic}]

Here are the six real valued roots to the system:

roots