Iterative Solvers: 1D aka ℝ

We want to solve non-linear equations of the form  y=(Subscript[y, 0])+h f[y] and other similar equations with small h.  We would like the solver to be general, fast, robust, reliable and cheap.  The problem is that a solver can not really be all of the above!

The first thing to know is if there is a solution: 
1) Clearly there is a solution with h=0 it is y=(Subscript[y, 0]).
2) If I pretend I have solved for y as a function of h i.e. y[h] I get
	y[h]=(Subscript[y, 0])+h f[y[h]]
and I can differentiate w.r.t h to get.
	y'[h]=f[y[h]]+h f'[y[h]]y'[h]
which I can solve to get 
	y'[h]=f[y[h]]/(1-h f'[y[h]])
	y[0]==(Subscript[y, 0])
provided the denominator is non-zero.  This is an ODE we know has a solution for small h because as h->0 the denominator goes to 1..
Conclusion:  There is a solution for small h  under reasonable smoothness conditions.

The simplest way to try to find y[h] for a small h  is a fixed point iteration
	(Subscript[y, k+1])=(Subscript[y, 0])+h f[(Subscript[y, k])]
if it converges the limit clearly satisfies our equation.

Newton's method should also work.  It is usaually described as solving g[y]==0 with the iteration
	(Subscript[y, k+1])=(Subscript[y, k])-g[y]/g'[y]
which translates (with g[y]=y-h f[y]-y0) into
	(Subscript[y, k+1])=(Subscript[y, k])+((Subscript[y, k])-h f[(Subscript[y, k])]-y0)/(1- h f'[(Subscript[y, k])])
If it converges this should converge faster than the fixed point iteration

Lazy Newton says skip evaluating the derivative when you can so
	(Subscript[y, k+1])=(Subscript[y, k])+((Subscript[y, k])-h f[(Subscript[y, k])]-y0)/(1- h f'[y0])
If it converges this should converge faster than the fixed point iteration

Example

f[y_]:= Sin[y](y^2)+ Cos[(y^2)];
y0=1.2;
CurvePic=ContourPlot[ y0+h f[y]- y,{y, 0, 3},{h,-1,1},
Contours->{0},
Epilog->{Red,PointSize[0.02],Point[{y0,0}]},
FrameLabel->{"y","h"}]

Fixed Point Iteration

y=y0;
hValFP=0.28;
FPPic=ListPlot[
FPData=Join[
{{y0,hValFP}},
Table[
y = y0+hValFP f[y];
{y,hValFP},
{1000}]
],
Joined->True,
PlotRange->All,
Mesh->All,
PlotStyle->Green
]

ListPlot[FPData[[All,1]],
PlotRange->All]

Newton

y=y0;
hValN=0.47;
NewtPic=ListPlot[
NewtData=Join[
{{y0,hValN}},
Table[
y = y - (y-hValN f[y]-y0)/(1-hValN f'[y]);
{y,hValN},
{10}]
],
Joined->True,
PlotRange->All,
Mesh->All,
PlotStyle->Cyan
]

ListPlot[NewtData[[All,1]],
PlotRange->All]

Lazy Newton

y=y0;
hValLazyN=0.48;
derval=f[y0]
LazyNewtPic=ListPlot[
LazyNewtData=Join[
{{y0,hValLazyN}},
Table[
y = y - (y-hValLazyN f[y]-y0)/(1-hValLazyN derVal);
{y,hValLazyN},
{i,0,100}]
],
Joined->True,
PlotRange->All,
Mesh->All,
PlotStyle->Pink
]

ListPlot[LazyNewtData[[All,1]],
PlotRange->All]

Combining Pics

Show[CurvePic,FPPic,LazyNewtPic, NewtPic]

Conclusions for 1D test cases.

Iterative Solvers: (ℝ^n)

All that changes is that I put vectors on top of everything and I need to solve a linear system instead of dividing.   
FP: 	(Subscript[y, k+1])=(Subscript[y, 0])+h f[(Subscript[y, k])]
N:	(Subscript[y, k+1])=(Subscript[y, k])+LinearSolve[Id - h ∇f[(Subscript[y, k])],((Subscript[y, k])-h f[(Subscript[y, k])]-y0)]
LN:	(Subscript[y, k+1])=(Subscript[y, k])+LinearSolve[Id - h ∇f[y0],((Subscript[y, k])-h f[(Subscript[y, k])]-y0)]
and we could of course update the frozen gradient value in LN if/when we wanted to.

Example 2D

f[y_]:= {
Sin[y](y^2)+ Cos[(y^2)],
y Log[ 1+ (y^2)]
}
y0={1.2,1.0};
y=y0;
hvalFP=0.1;
Table[
y = y0+hValFP f[y],
{100}]

Fixed Point Iteration

y=y0;
hValFP=0.28;
FPPic=ListPlot[
FPData=Join[
{{y0,hValFP}},
Table[
y = y0+hValFP f[y];
{y,hValFP},
{1000}]
],
Joined->True,
PlotRange->All,
Mesh->All,
PlotStyle->Green
]

ListPlot[FPData[[All,1]],
PlotRange->All]

Newton

y=y0;
hValN=0.47;
NewtPic=ListPlot[
NewtData=Join[
{{y0,hValN}},
Table[
y = y - (y-hValN f[y]-y0)/(1-hValN f'[y]);
{y,hValN},
{10}]
],
Joined->True,
PlotRange->All,
Mesh->All,
PlotStyle->Cyan
]

ListPlot[NewtData[[All,1]],
PlotRange->All]

Lazy Newton

y=y0;
hValLazyN=0.48;
derval=f[y0]
LazyNewtPic=ListPlot[
LazyNewtData=Join[
{{y0,hValLazyN}},
Table[
y = y - (y-hValLazyN f[y]-y0)/(1-hValLazyN derVal);
{y,hValLazyN},
{i,0,100}]
],
Joined->True,
PlotRange->All,
Mesh->All,
PlotStyle->Pink
]

ListPlot[LazyNewtData[[All,1]],
PlotRange->All]

Combining Pics

Show[CurvePic,FPPic,LazyNewtPic, NewtPic]

Conclusions for 1D test cases.