Quick guide to Mathematica

To get help, place the cursor on the function and press F1
To run all commands in one cell, place the cursor there and press Shift + Enter or Enter on the numerical keyboard.
Output of many commands is suppressed by ; at the end of the command. If you want to see the output, delete ;
Argument of functions must be in square brackets [...], braces {...} are for arrays, ranges etc., use double brackets [[...]] to access elements of arrays
expr /. {x -> a} means substitute a for x in expr, expr // function means apply function to expr, i.e. it is equivalent to function[expr].
D[f[x],x] = derivative of f[x] with respect to x.
To define a function of x one can use f[x_] := 1+(x^2), notice the underscore and the colon, or # and & like in f := 1+((#)^2) &
N[expr] evaluates expr with machine precision, if you use the decimal point in the expression it will be also evaluated with machine precision
Special characters can be inserted by pressing Esc ... Esc, e.g. Esc p Esc gives π
Clear[...] is used to unset any variables which could have been assigned previously.
Some other useful commands: Simplify, Expand, Factor; Integrate, Series, Sum

Clear all symbols from previous evaluations to avoid problems

Clear["Global`*"]

Solution of algebraic equations

Basic iterations

To solve an equation

x=f(x)

iteratively we write one iteration step as

(Subscript[y, k+1])=f((Subscript[y, k]))

Let us apply this simple method to an equation

x=cos(x-α)

and study convergence with respect to the parameter α ∈[0,1].

α=0;
f[z_]:=Cos[z-α];
n=30;
x0=0.0;
(* simple iterations *)
x=ConstantArray[0.0,n];x[[1]]=x0;
Do[x[[i]]=f[x[[i-1]]],{i,2,n}];
Print["Simple iterations: ",x]
(* Mathematica's numerical solution *)
root=t/.FindRoot[t==f[t],{t,x0}];
Print["Mathematica's solution: ", root]
(* plot relative error of the consequent approximations *)
ListLogPlot[{Abs[root-x]/root},AxesLabel->{"n","Relative error"},PlotLegends->{"Relative error of simple iterations"}]

plotFunc=Plot[{f[x],x},{x,-0.2,1.2},AspectRatio->1,AxesLabel->{"x","f(x)"}];plotIter=Graphics[Table[Line[{{x[[i]],x[[i]]},{x[[i]],x[[i+1]]},{x[[i+1]],x[[i+1]]}}],{i,1,n-1}]];
Show[plotFunc,plotIter]

Acceleration of convergence - Aitken’s (Δ^2) process

For sequences which converge linearly we can accelerate convergence by creating a new sequence by

(Subscript[y, n])=(Subscript[x, n])-(((((Subscript[Δx, n])))^2)/((Δ^2)(Subscript[x, n]))),

where

(Subscript[Δx, n]) = (Subscript[x, n+1])-(Subscript[x, n])  and  (Δ^2)(Subscript[x, n]) = ((Subscript[x, n+2])-(Subscript[x, n+1]))-((Subscript[x, n+1])-(Subscript[x, n]))=(Subscript[x, n+2])-2(Subscript[x, n+1])+(Subscript[x, n])

are forward differences. We can repeat this process again and again.
These new sequences converge faster but, in general, also linearly.

Ax1=MyAitkenDeltaSquareProcess[x];
Ax2=MyAitkenDeltaSquareProcess[Ax1];
Ax3=MyAitkenDeltaSquareProcess[Ax2];
Ax4=MyAitkenDeltaSquareProcess[Ax3];
ListLogPlot[{Abs[root-x],Abs[root-Ax1],Abs[root-Ax2],Abs[root-Ax3],Abs[root-Ax4]},PlotStyle->{Blue,Green,Red,Brown,Black},AxesLabel->{"n","Relative error"},PlotLegends->{"Simple iterations","First Aitken (Δ^2) accelaration","Second Aitken (Δ^2) accelaration","Third Aitken (Δ^2) accelaration","Fourth Aitken (Δ^2) accelaration"}]

Acceleration of convergence - Aitken-Steffensen method

A method for acceleration of linear convergence of simple iterations when solving x=f(x) by introducing a new series

(Subscript[y, k])=(Subscript[x, k])-((((f( (Subscript[x, k]))-(Subscript[x, k])))^2)/(f(f((Subscript[x, k])))-2f((Subscript[x, k]))+(Subscript[x, k])))

This new sequences converge quadratically .

Sx=MyAitkenSteffensenFindRoot[f,x0,n];
ListLogPlot[{Abs[root-x],Abs[root-Ax1],Abs[root-Ax2],Abs[root-Ax3],Abs[root-Sx]},PlotStyle->{Blue,Green,Green,Green,Red},PlotRange->{10^-20,1},PlotRange->{10^-20,1},AxesLabel->{"n","Relative error"},PlotLegends->{"Simple iterations","First to third Aitken (Δ^2) accelaration",None,None,"Aitken-Steffensen method"}]

Methods to find roots

To find a root of a function g(x), i.e. a solution of

g(x)=0,

we can use some standard method, such as the bisection or regula falsi method, or we can write

x=f(x)=x-ϕ(x)g(x)

and find the root iteratively as above. 
By choosing  ϕ(x) as 1/ϕ’(x), we get Newton’s method, which converges quadratically.
By approximating the derivative with the difference, we get secant method which converges also quadratically.
The last compared method is consequent Pade [1/1] approximations.

Comparison of bisection, Newton’s, secant, and Pade’s methods

The first example, solution of x = cos(x):

g[z_]=z-Cos[z];
Dg[z_]=D[g[z],z];
n=20;a=0.0; b=1.0; c=0.7;
xB=MyBisectionFindRoot[g,a,b,n];
xN=MyNewtonFindRoot[g,Dg,a,n];
xS=MySecantFindRoot[g,a,b,n];
xP=MyPadeFindRoot[g,a,b, c,n];
root=t/.FindRoot[t==Cos[t],{t,a}];
ListLogPlot[{Abs[root-xB],Abs[root-xN],Abs[root-xS],Abs[root-xP]},PlotStyle->{Blue,Green,Red, Black},
AxesLabel->{"n","Relative error"},PlotLegends->{"Bisection method","Newton's method","Secant method","Pade's method"}]

The second example, solution of (x^p) = a, a real root:

p=3;
a=3;
g[z_]=z^p-3;
Dg[z_]=D[g[z],z];
n=20;a=1.3; b=2.0; c=1.7;
xB=MyBisectionFindRoot[g,a,b,n];
xN=MyNewtonFindRoot[g,Dg,a,n];
xS=MySecantFindRoot[g,a,b,n];
xP=MyPadeFindRoot[g,a,b, c,n];
root=3^(1/p);
ListLogPlot[{Abs[root-xB],Abs[root-xN],Abs[root-xS],Abs[root-xP]},PlotStyle->{Blue,Green,Red, Black},AxesLabel->{"n","Relative error"},PlotLegends->{"Bisection method","Newton's method","Secant method","Pade's method"}]

The second example, solution of (x^p) = a, a complex root:

p=3;
a=3;
g[z_]=z^p-3;
Dg[z_]=D[g[z],z];
n=20;a=-1.0-I; b=0.0; c = -I;
xN=MyNewtonFindRoot[g,Dg,a,n];
xS=MySecantFindRoot[g,a,b,n];
xP=MyPadeFindRoot[g,a,b, c,n];
root=3^(1/3)Exp[-2I π/3];
ListLogPlot[{Abs[root-xN],Abs[root-xS],Abs[root-xP]},PlotStyle->{Blue,Green,Red, Black},AxesLabel->{"n","Relative error"},PlotLegends->{"Newton's method","Secant method","Pade's method"}]

Newton’s fractal

When solving (z^3) = 1 by Newton’s method, we can end up at different complex roots depending where we start.
Moreover, the number of iterations to achieve a given tolerance also depends significantly on the starting point

g[z_]=(z^3)-1;
Dg[z_]=D[g[z],z];
tol=1.0*^-8;
itermax=100;
Do[
z0=-0.01+I*(0.0050699+i*0.00000001);
{zr,niter}=MyNewtonFindRootWithConvergenceTest[g,Dg,z0,tol,itermax];
Print["Root ",zr," found after ",niter," iterations"],
{i,1,9}
];

Newton’s fractal is obtained when each point in the complex plane is coloured according to which root of (z^3) = 1 Newton’s method converges when starting at that point.

{nx,ny}={500,500};
{xmin,xmax}={-0.2,0.2};
{ymin,ymax}={-0.2,0.2};
dx = (xmax-xmin)/(nx-1);
dy = (ymax-ymin)/(ny-1);
X = Range[xmin, xmax, dx];
Y = Range[ymin, ymax, dy];
colors=ConstantArray[0,{ny,nx}];
tol=10.0*^-6;
nmax=50;
Do[
Do[
a=X[[i]]+I Y[[j]];
{xr,nc}=MyNewtonFindRootWithConvergenceTest[g,Dg,a,tol,nmax];
Which[Abs[xr-1]<tol,colors[[ny-j+1,i]]=nc,
Abs[xr-(-1/2+I Sqrt[3]/2)]<tol,colors[[ny-j+1,i]]=100+nc,
Abs[xr-(-1/2-I Sqrt[3]/2)]<tol,colors[[ny-j+1,i]]=200+nc,
True,colors[[ny-j+1,i]]=0
],
{j,1,ny}
],
{i,1,nx}
];
MyColor[ a_ ]:=Which[
a==0,Black,
a<100,RGBColor[(100-3(a-10))/100,0,0],
a<200,RGBColor[0,(100-3(a-110))/100,0],
a≥200,RGBColor[0,0,(100-3(a-210))/100]];
ArrayPlot[colors,ColorFunctionScaling->False,
ColorFunction->MyColor]

x = (-5, 5), y = (-5, 5)                                                              x = (0.3, 0.8), y = (1.0, 1.5)

x = (0.4, 0.45)                                    x = (0.4, 0.41)                                  x = (0.4, 0.401)
y = (1.2, 1.25)                                    y = (1.2, 1.21)                                  y = (1.2005, 1.2015)

Finding roots of vector functions of more variables

In general, we can use Newton’s method to find roots of g(x) where g is a vector function of n variables with n components.
Here is an example for two functions of two variables.

Clear[x,y]
f1[x_,y_]=Exp[x-2]-y;
f2[x_,y_]=-x+y^2;
Plot3D[{f1[x,y],f2[x,y],0},{x,-1,1},{y,-1,1}]
rootMath={x,y}/.FindRoot[{f1[x,y]==0,f2[x,y]==0},{{x,1},{y,1}}];
Print["Mathematica's solution: ", rootMath]

Iterations by Newton’s method and comparison with Mathematica’s solution

f[{x_,y_}]={f1[x,y],f2[x,y]};
J[{x_,y_}]={{D[f1[x,y],x],D[f1[x,y],y]},{D[f2[x,y],x],D[f2[x,y],y]}};
Print["Jacobian = ",J[{x,y}]//MatrixForm]
x0={0.5,1.0};
rootNewtonND=MyNewtonFindRootND[f,J,x0,10];
ListLogPlot[Map[Norm[#-rootMath]&,rootNewtonND],PlotStyle->{Blue},AxesLabel->{"n","Absolute error"},PlotLegends->{"Newton's method"}]
error=Map[Norm[#-rootMath]&,rootNewtonND];

Checking the result

Map[f,xNewtonND]