{NotebookFileName[],DateString[]}

style[
t_
]:=
Style[t,Bold,FontFamily->"Arial",FontSize->14]

newtonStep

newtonStep[f] is a function that, applied to an estimate (Subscript[x, i]), will take a Newton’s method step (towards, presumably, a root of f). It can be used as the first argument to Nest or NestList to do a Newton iteration.

newtonStep[
f_
]:=
Function[xi,
Block[{m,x,y},
m=D[f[x],x]/.x->xi;
y=f[xi];
xi-(y/m)
]
]

Same as newtonStep, but guaranteed to return a machine precision number:

newtonStepN[
f_
]:=
Function[xi,
Block[{m,x,y},
m=D[f[x],x]/.x->xi;
y=f[xi];
N[xi-(y/m)]
]
]

newtonStep[x|->(x^2)-2]

newtonStep[x|->(x^2)-2][2]

NestList[newtonStep[x|->(x^2)-2],2,4]

NestList[newtonStepN[x|->(x^2)-2],2,4]

NestList[newtonStep[x|->(x^3)],1,4]

NestList[newtonStep[x|->(x^p)],1,4]//Simplify

Solve[((p-1)/p)==1]

Plot[((p-1)/p),{p,-10,10}
]

bisectionStep

bisectionStep[
f_
]:=
Function[i,
Block[{xl,xr,xm,fl,fr,fm},
{xl,xr}=i;
xm=((xl+xr)/2);
{fl,fr,fm}=f/@{xl,xr,xm};
If[fl*fr>0,
Abort[];
];
If[fl*fm<0,
{xl,xm},
{xm,xr}
]
]
]

bisectionStepN[
f_
]:=
Function[i,
Block[{xl,xr,xm,fl,fr,fm},
{xl,xr}=i;
xm=N[((xl+xr)/2)];
{fl,fr,fm}=N[f/@{xl,xr,xm}];
If[fl*fr>0,
Abort[];
];
If[fl*fm<0,
{xl,xm},
{xm,xr}
]
]
]

bisectionStep[x|->(x^2)-2]

bisectionStep[x|->(x^2)-2][{1,2}]

NestList[bisectionStep[x|->(x^2)-2],{1,2},4]

bisectionStepN[x|->(x^2)-2][{1,2}]

NestList[bisectionStepN[x|->(x^2)-2],{1,2},4]

secantStep

secantStep[
f_
]:=
Function[last2,
Block[{x1,x0,f1,f0x},
{x0,x1}=last2;
{f0,f1}=f/@{x0,x1};
x=x1-f1(((x1-x0)/(f1-f0)));
{x1,x}
]
]

secantStepN[
f_
]:=
Function[last2,
Block[{x1,x0,f1,f0,x},
{x0,x1}=N[last2];
{f0,f1}=N[f/@{x0,x1}];
x=N[x1-f1(((x1-x0)/(f1-f0)))];
{x1,x}
]
]

secantStep[x|->(x^2)-2]

secantStep[x|->(x^2)-2][{1,2}]

NestList[secantStep[x|->(x^2)-2],{1,2},4]

secantStepN[x|->(x^2)-2][{1,2}]

NestList[secantStepN[x|->(x^2)-2],{1,2},4]

Pendulum

$Assumptions={{θ,θ0}∈Reals}

{(θ|θ0)∈Reals,0<θ<π,0<θ0<π}

Integrate[(1/(Sqrt[Cos[θ]-Cos[θ0]])),{θ,0,θ0},GenerateConditions->False]

tθ[
θ0_,
L_,
g_
]:=Evaluate[
4(Sqrt[(L/(2g))])Integrate[(1/(Sqrt[Cos[θ]-Cos[θ0]])),{θ,0,θ0},GenerateConditions->False]
];
tθ[θ0,L,g]

Limit[tθ[θ0,L,g],θ0->0]

Series[tθ[θ0,1,1],{θ0,0,3}]

Series[tθ[θ0,1,1],{θ0,θ1,3}]

D[EllipticF[((θ0)/2),((Csc[((θ0)/2)])^2)],θ0]

D[EllipticF[((θ0)/2),Cθ02],θ0]

D[EllipticF[θ02,((Csc[((θ0)/2)])^2)],θ0]

Old example

f2[
x_
]:=
(E^x)+Log[x]

NestList[bisectionStep[f2],{0,1},2]

NestList[bisectionStepN[f2],{0,1},5]

Fortunately for me, Mathematica is able to evaluate log 0 and compare it, so my bisection step doesn’t abort on [0, 1].

Log[0]*E>0

NestList[newtonStep[f2],0.5,5]

ϕ==GoldenRatio/.Solve[1+ϕ==(ϕ^2),ϕ]

NestList[newtonStep[x|->Log[x]],(1/2),4]

NestList[newtonStep[x|->Log[x]],0.5`20,4]

Mobius lecture examples

fr1[
x_
]:=
(x^5)-3

{{1,2},fr1/@{1,2}}

x=.

frline1=ListInterpolation[fr1/@{1,2},{1,2},InterpolationOrder->1]

frint1=x/.N[Solve[frline1[x]==0,x][[1]]]

Plot[frline1[x],{x,0.8,2.2},
Epilog->{
PointSize[Large],
Point[{{1,fr1[1]},{2,fr1[2]}}],
Red,
Point[{{frint1,frline1[frint1]}}]
}
]

Block[{xmin=1.08,xmax=1.3,x1=1.1,m1,p1,x2,p2},
m1=fr1'[x1];
p1={x1,fr1[x1]};
x2=newtonStep[fr1][x1];
p2={x2,0};
Plot[fr1[x],{x,xmin,xmax},
Epilog->{
PointSize[Large],
Point[{p1,p2}],
Line[{{p1,p2}}],
Text[style["(1.1, f(1.1))"],p1,{-1.3,0}]
}
]
]

Block[{xmin=1.08,xmax=1.3,x1=1.1,m1,p1,x2,p2,p2f,x3,p3},
m1=fr1'[x1];
p1={x1,fr1[x1]};
x2=newtonStep[fr1][x1];
p2={x2,0};
p2f={x2,fr1[x2]};
x3=newtonStep[fr1][x2];
p3={x3,0};
Plot[fr1[x],{x,xmin,xmax},
Epilog->{
PointSize[Large],
Point[{p1,p2,p2f,p3}],
Red,
Line[{p2,p2f}],
Black,
Line[{{p1,p2}}],Line[{p2f,p3}],
Text[style["(1.1, f(1.1))"],p1,{-1.3,0}],
Text[style["(Subscript[x, 2])"],p2,{1,1}]
}
]
]

x1=1.1`20

fr1[x1]

fr1'[x1]

fr1sols=NestList[newtonStep[fr1],x1,4]

((3.0`20)^(1/5))

fr1/@fr1sols

fr1sols-(3^(1/5))

((fr1sols(fr1/@fr1sols))/5)

((fr1/@fr1sols)/3)

((fr1/@fr1sols)/(3 5))

((fr1/@fr1sols)/(3 5))fr1sols

(fr1sols^5)

Why are these not exactly equal? Let’s take (Subscript[x, 1])=1.1, where the difference is biggest.

f((Subscript[x, 1]))=(Subscript[x, 1]^5)-3
(Subscript[x, 2])=(Subscript[x, 1])-((f((Subscript[x, 1])))/(f'((Subscript[x, 1]))))

(((fr1/@fr1sols)/(3 5))fr1sols)((3/(fr1sols^5)))

Differences[fr1sols]

OK, that’s it. They’re not equal because the (Subscript[x, i]^5) are not exactly equal to 3.

((fr1/@fr1sols)/(fr1'/@fr1sols))

Yes, those are exactly the differences.

((f((Subscript[x, 1])))/(f'((Subscript[x, 1]))))=(((Subscript[x, 1]^5)-3)/(5(Subscript[x, 1]^4)))=((Subscript[x, 1])/5)-((3(Subscript[x, 1]))/(Subscript[x, 1]^5))

fr1'/@fr1sols

x2=(3^(1/5));
x3=x2+(10^(-5))

((fr1[x3])/(3 5))//N

(((x3^5)-(x2^5))/(x2^5))//N

((x3-x2)/x2)//N

((((x3^5)-(x2^5))/(x2^5))/((x3-x2)/x2))//N

fr1sol2=NestList[bisectionStep[fr1],{1.0`20,1.4`20},10]

Mean/@fr1sol2

Mean[fr1sol2[[-1]]]

fr1[Mean[fr1sol2[[-1]]]]

fm1[
x_
]:=
Sin[x]-(E^x)+2

fm1'[x]

1-((fm1[1])/(fm1'[1]))//N

newtonStep[fm1]

NestList[newtonStepN[fm1],1,4]

fm2[
x_
]:=
Cos[x]-x

NestList[bisectionStepN[fm2],{0,1},5]

NestList[newtonStepN[fm2],0.5,5]

Sum of exponentials

I’m thinking of this as a HIF degradation model. I have (0^th)-order production of HIF, followed by first-order proline oxidation, followed by first-order degradation of  oxidized HIF. I just pulled the kinetic constants out of my ass.

Hprod=1;
kOx=1;
kDeg=0.115;
O2=.

eqsH={
H'[t]== Hprod-kOx O2 H[t],
Ho'[t]==kOx H[t]-kDeg Ho[t]
}

icsH={
H[0]==H0,Ho[0]== Ho0
}

solH=DSolve[Join[eqsH,icsH],{H[t],Ho[t]},t][[1]]//Simplify

H[t]+Ho[t]/.solH

H[t]+Ho[t]/.solH/.O2->0.1//Simplify

Solve[0.5 Hmax==H[t]+Ho[t]/.solH,t,Reals]

Solve[0.9 Hmax==H[t]+Ho[t]/.solH,t,Reals]

H[t]+Ho[t]/.solH/.O2->0.21

Hs0=Thread[{H0,Ho0}->Limit[{H[t],Ho[t]}/.solH/.O2->0.21,t->∞]]

hypO2=0.05

H[t]+Ho[t]/.solH/.Hs0/.O2->hypO2//Simplify

Plot[Evaluate[
{Ho[t],H[t],H[t]+Ho[t]}/.solH/.Hs0/.O2->hypO2
],{t,0,100}
,PlotLegends->{(Subscript[H, o])[t],H[t],H[t]+(Subscript[H, o])[t]}
]

Hhyp[
t_
]:=Evaluate[
H[t]+Ho[t]/.solH/.Hs0/.O2->hypO2//Simplify
];
Hhyp[t]

Clear[Hhypi];
Hhypi[
t_
]:=
194+102(E^(-0.115`20t))-250(E^(-0.05`20t));
Hhypi[t]

Hsimp[
t_
]:=Evaluate[
Hhyp[t][[{1,3}]]
];
Hsimp[t]

Hsimpi[
t_
]:=Evaluate[
Hhypi[t][[{1,3}]]
];
Hsimpi[t]

97/250//N

Log[97/250//N]

((Log[97/250//N])/(-0.05))

Hhypi[0]

Plot[Hhypi[t],{t,0,100}
,PlotRange->{All,{0,All}}
,PlotLegends->{"[H]"[t]}
,AxesLabel->{style["t (min)"],None}
]

Hmax=Limit[Hhypi[t],t->∞]

Solve[Hsimpi[t]==0.5Hmax,t,Reals]

t1=SetPrecision[t/.Solve[Hsimpi[t]==0.5Hmax,t,Reals][[1]],20]

Solve[Hhyp[t]==0.5Hmax,t,Reals]

Solve[Hsimp[t]==0.9Hmax,t,Reals]

Solve[Hhyp[t]==0.9Hmax,t,Reals]

f3[
t_
]:=
Hhypi[t]-(1/2)Hmax;
f3[t]

f3'[t]

{Hhypi[t1],f3[t1]}

f3sols=NestList[newtonStep[f3],t1,6]

Differences[f3sols]

f3'/@f3sols

f3/@f3sols

Block[{tmin=15,tmax=19,m1,p1,t2,p2,p2f,t3,p3},
m1=f3'[t1];
p1={t1,f3[t1]};
t2=newtonStep[f3][t1];
p2={t2,0};
p2f={t2,f3[t2]};
t3=newtonStep[f3][t2];
p3={t3,0};
Plot[f3[t],{t,tmin,tmax}
,AxesLabel->{style["t (min)"],None}
,Epilog->{
PointSize[Large],
Point[{p1,p2}],
Line[{{p1,p2}}],
Text[style["((Subscript[t, 1]), f((Subscript[t, 1]))"],p1,{1.3,0}],
Text[style["(Subscript[t, 2])"],p2,{1,-1}]
}
]
]

Block[{tmin=15,tmax=16,m1,p1,t2,p2,p2f,t3,p3},
m1=f3'[t1];
p1={t1,f3[t1]};
t2=newtonStep[f3][t1];
p2={t2,0};
p2f={t2,f3[t2]};
t3=newtonStep[f3][t2];
p3={t3,0};
Plot[f3[t],{t,tmin,tmax}
,AxesLabel->{style["t (min)"],None}
,Epilog->{
PointSize[Large],
Point[{p1,p2,p2f,p3}],
Red,
Line[{p2,p2f}],
Black,
Line[{{p1,p2}}],Line[{p2f,p3}],
Text[style["((Subscript[t, 1]), f((Subscript[t, 1]))"],p1,{1.3,0}],
Text[style["(Subscript[t, 2])"],p2,{1,-1}],
Text[style["(Subscript[t, 3])"],p3,{1,-1}]
}
]
]

Plot[Hhypi[t],{t,0,100}
,AxesLabel->{style["t (min)"], None}
,PlotRange->{All,{0,All}}
,PlotLegends->{"[H]"[t]}
,Epilog->{
PointSize[Large],
Point[{{f3sols[[5]],Hhypi[f3sols[[5]]]}}],
Text[style["((Subscript[t, 5]), 97)"],{f3sols[[5]],Hhypi[f3sols[[5]]]},{-1.4,0}]
}
]

t1

f3[t1]

t2=t1-1

Block[{tmin=15,tmax=19,m1,p1,p2,p2f,t3,p3},
p1={t1,f3[t1]};
p2={t2,f3[t2]};
m=((f4[t2]-f3[t1])/(t2-t1));
t3=secantStep[f3][{t2,t1}][[2]];
p3={t3,0};
Plot[f3[t],{t,tmin,tmax}
,AxesLabel->{style["t (min)"], None}
,Epilog->{
PointSize[Large],
Point[{p1,p2,p3}],
Line[{{p1,p3}}],
Text[style["((Subscript[t, 1]), f((Subscript[t, 1])))"],p1,{1.3,0}],
Text[style["((Subscript[t, 2]), f((Subscript[t, 2])))"],p2,{1.3,0}],
Text[style["(Subscript[t, 3])"],p3,{1.3,-1}]
}
]
]

f3ssols=NestList[secantStep[f3],{t2,t1},8]

Last/@f3ssols

f3sols

Limitations

x=.

f4[
x0_:0
]:=
With[{d=1+x0},
Function[x,
x (E^x)-(E^x)+d
]
];
f4[][x]

f4[][0]

Plot[Evaluate[f4[0.01][x]],{x,-10,3}
,PlotLegends->{x (E^x)-(E^x)+1.01}
]

NestList[newtonStep[f4[0.01]],0.5,5]

Block[{is=400,xmin=-1,xmax=1.2,x0=0.01,x1=-0.5,nsteps=5,steps,curvepoints,axispoints,labels,colors},
steps=NestList[newtonStep[f4[x0]],x1,nsteps];
curvepoints=Transpose[{steps,f4[x0]/@steps}];
axispoints=Transpose[{steps,ConstantArray[0.0,nsteps+1]}];
colors=Table[ColorData[97][i],{i,nsteps+1}];
labels=Table[
Text[style[ToString[i]],curvepoints[[i]],{0,-1}],
{i,nsteps+1}
];
labels=Riffle[colors,labels];
axispoints=Riffle[colors,Point[{#}]&/@axispoints];
curvepoints=Riffle[colors,Point[{#}]&/@curvepoints];
Plot[f4[x0][x],{x,xmin,xmax}
,PlotRange->{All,All}
,PlotLegends->{f4[x0][x]}
,ImageSize->is
,Epilog->{
PointSize[Large],
curvepoints,
axispoints,
labels
}
]
]

Block[{is=400,xmin=-1,xmax=1.2,x0=0.01,x1=-0.5,nsteps=5,steps,curvepoints,axispoints,labels,colors},
steps=NestList[newtonStep[f4[x0]],x1,nsteps];
curvepoints=Transpose[{steps,f4[x0]/@steps}];
axispoints=Transpose[{steps,ConstantArray[0.0,nsteps+1]}];
colors=Table[ColorData[97][i],{i,nsteps+1}];
labels=Table[
Text[style[ToString[i]],curvepoints[[i]],{0,-1}],
{i,nsteps+1}
];
labels=Riffle[colors,labels];
axispoints=Riffle[colors,Point[{#}]&/@axispoints];
curvepoints=Riffle[colors,Point[{#}]&/@curvepoints];
Plot[f4[x0][x],{x,xmin,xmax}
,PlotRange->{All,{-0.01,0.1}}
,PlotLegends->{f4[x0][x]}
,ImageSize->is
,Epilog->{
PointSize[Large],
curvepoints,
axispoints,
labels
}
]
]

Plot[Evaluate[f4[0.01][x]],{x,-10,3}
,PlotLegends->{x (E^x)-(E^x)+1}
]

Block[{is=400,xmin=-0.5,xmax=0.1,x0=0,x1=-0.5,nsteps=5,steps,curvepoints,axispoints,labels,colors},
steps=NestList[newtonStep[f4[x0]],x1,nsteps];
curvepoints=Transpose[{steps,f4[x0]/@steps}];
axispoints=Transpose[{steps,ConstantArray[0.0,nsteps+1]}];
colors=Table[ColorData[97][i],{i,nsteps+1}];
labels=Table[
Text[style[ToString[i]],curvepoints[[i]],{0,-1}],
{i,nsteps+1}
];
labels=Riffle[colors,labels];
axispoints=Riffle[colors,Point[{#}]&/@axispoints];
curvepoints=Riffle[colors,Point[{#}]&/@curvepoints];
{
Plot[f4[x0][x],{x,xmin,xmax}
,PlotRange->{All,All}
,PlotLegends->{f4[x0][x]}
,ImageSize->is
,Epilog->{
PointSize[Large],
curvepoints,
axispoints,
labels
}
],
steps
}
]//Column

(((-0.5))^(1/3))

Surd[-0.5,3]

Solve[Surd[x,3]-0.9==0,x]

(0.9^3)

Block[{is=400,xmin=-0.5,xmax=0.1,x0=0,x1=-0.5,nsteps=5,f,steps,curvepoints,axispoints,labels,colors},
f=(x|->Surd[x,3]-0.5);
steps=NestList[newtonStep[f],x1,nsteps];
xmin=Min[Append[steps,xmin]];
xmax=Max[Append[steps,xmax]];
curvepoints=Transpose[{steps,f/@steps}];
axispoints=Transpose[{steps,ConstantArray[0.0,nsteps+1]}];
colors=Table[ColorData[97][i],{i,nsteps+1}];
labels=Table[
Text[style[ToString[i]],curvepoints[[i]],{0,-1}],
{i,nsteps+1}
];
labels=Riffle[colors,labels];
axispoints=Riffle[colors,Point[{#}]&/@axispoints];
curvepoints=Riffle[colors,Point[{#}]&/@curvepoints];
{
Plot[f[x],{x,xmin,xmax}
,PlotRange->{All,All}
,PlotLegends->{f[x]}
,ImageSize->is
,Epilog->{
PointSize[Large],
curvepoints,
axispoints,
labels
}
],
steps
}
]//Column

Block[{is=400,xmin=-0.3,xmax=0.1,x0=0,x1=-0.05,nsteps=5,f,steps,curvepoints,axispoints,labels,colors},
f=(x|->Surd[x,3]-0.5);
steps=NestList[newtonStep[f],x1,nsteps];
xmin=Min[Append[steps,xmin]];
xmax=Max[Append[steps,xmax]];
curvepoints=Transpose[{steps,f/@steps}];
axispoints=Transpose[{steps,ConstantArray[0.0,nsteps+1]}];
colors=Table[ColorData[97][i],{i,nsteps+1}];
labels=Table[
Text[style[ToString[i]],curvepoints[[i]],{0,-1}],
{i,nsteps+1}
];
labels=Riffle[colors,labels];
axispoints=Riffle[colors,Point[{#}]&/@axispoints];
curvepoints=Riffle[colors,Point[{#}]&/@curvepoints];
{
Plot[f[x],{x,xmin,xmax}
,PlotRange->{All,All}
,PlotLegends->{f[x]}
,ImageSize->is
,Epilog->{
PointSize[Large],
curvepoints,
axispoints,
labels
}
],
steps
}
]//Column

Plot[f4[-0.001][x],{x,-0.5,0.5}
,PlotRange->{-0.002,0.1}
]

Plot[Evaluate[f4[-0.001][x]],{x,-10,3}
,PlotLegends->{f4[-0.001][x]}
]

Block[{x0=-0.001,roots},
roots={
x/.FindRoot[f4[x0][x],{x,0.05}],
x/.FindRoot[f4[x0][x],{x,-0.05}]
};
roots
]

Block[{is=400,xmin=-0.5,xmax=1.0.5,x0=-0.001,x1=0.5,nsteps=5,steps,curvepoints,axispoints,labels,colors,roots,rootpoints},
roots={
x/.FindRoot[f4[x0][x],{x,0.05}],
x/.FindRoot[f4[x0][x],{x,-0.05}]
};
steps=NestList[newtonStep[f4[x0]],x1,nsteps];
xmin=Min[Append[steps,xmin]];
xmax=Max[Append[steps,xmax]];
curvepoints=Transpose[{steps,f4[x0]/@steps}];
axispoints=Transpose[{steps,ConstantArray[0.0,nsteps+1]}];
colors=Table[ColorData[97][i],{i,nsteps+1}];
labels=Table[
Text[style[ToString[i]],curvepoints[[i]],{0,-1}],
{i,nsteps+1}
];
labels=Riffle[colors,labels];
axispoints=Riffle[colors,Point[{#}]&/@axispoints];
curvepoints=Riffle[colors,Point[{#}]&/@curvepoints];
rootpoints=Point[Transpose[{roots,{0,0}}]];
{
Plot[f4[x0][x],{x,xmin,xmax}
,PlotRange->{All,{-0.01,All}}
,PlotLegends->{f4[x0][x]}
,ImageSize->is
,Epilog->{
PointSize[Large],
curvepoints,
axispoints,
labels,
{Black,rootpoints}
}
],
steps
}
]//Column

Block[{is=400,xmin=-0.5,xmax=1.0.5,x0=-0.001,x1=-0.5,nsteps=5,steps,curvepoints,axispoints,labels,colors,roots,rootpoints},
roots={
x/.FindRoot[f4[x0][x],{x,0.05}],
x/.FindRoot[f4[x0][x],{x,-0.05}]
};
steps=NestList[newtonStep[f4[x0]],x1,nsteps];
xmin=Min[Append[steps,xmin]];
xmax=Max[Append[steps,xmax]];
curvepoints=Transpose[{steps,f4[x0]/@steps}];
axispoints=Transpose[{steps,ConstantArray[0.0,nsteps+1]}];
colors=Table[ColorData[97][i],{i,nsteps+1}];
labels=Table[
Text[style[ToString[i]],curvepoints[[i]],{0,-1}],
{i,nsteps+1}
];
labels=Riffle[colors,labels];
axispoints=Riffle[colors,Point[{#}]&/@axispoints];
curvepoints=Riffle[colors,Point[{#}]&/@curvepoints];
rootpoints=Point[Transpose[{roots,{0,0}}]];
{
Plot[f4[x0][x],{x,xmin,xmax}
,PlotRange->{All,{-0.01,All}}
,PlotLegends->{f4[x0][x]}
,ImageSize->is
,Epilog->{
PointSize[Large],
curvepoints,
axispoints,
labels,
{Black,rootpoints}
}
],
steps
}
]//Column

Block[{is=400,xmin=-0.5,xmax=1.0.5,x0=-0.001,x1=-2,nsteps=5,steps,curvepoints,axispoints,labels,colors,roots,rootpoints},
roots={
x/.FindRoot[f4[x0][x],{x,0.05}],
x/.FindRoot[f4[x0][x],{x,-0.05}]
};
steps=NestList[newtonStep[f4[x0]],x1,nsteps];
xmin=Min[Append[steps,xmin]];
xmax=Max[Append[steps,xmax]];
curvepoints=Transpose[{steps,f4[x0]/@steps}];
axispoints=Transpose[{steps,ConstantArray[0.0,nsteps+1]}];
colors=Table[ColorData[97][i],{i,nsteps+1}];
labels=Table[
Text[style[ToString[i]],curvepoints[[i]],{0,-1}],
{i,nsteps+1}
];
labels=Riffle[colors,labels];
axispoints=Riffle[colors,Point[{#}]&/@axispoints];
curvepoints=Riffle[colors,Point[{#}]&/@curvepoints];
rootpoints=Point[Transpose[{roots,{0,0}}]];
{
Plot[f4[x0][x],{x,xmin,xmax}
,PlotRange->{All,{-0.01,All}}
,PlotLegends->{f4[x0][x]}
,ImageSize->is
,Epilog->{
PointSize[Large],
curvepoints,
axispoints,
labels,
{Black,rootpoints}
}
],
steps
}
]//Column

Block[{is=400,xmin=-0.5,xmax=1.0.5,x0=-0.001,x1=-8,nsteps=20,steps,curvepoints,axispoints,labels,colors,roots,rootpoints},
roots={
x/.FindRoot[f4[x0][x],{x,0.05}],
x/.FindRoot[f4[x0][x],{x,-0.05}]
};
steps=NestList[newtonStep[f4[x0]],x1,nsteps];
xmin=Min[Append[steps,xmin]];
xmax=Max[Append[steps,xmax]];
curvepoints=Transpose[{steps,f4[x0]/@steps}];
axispoints=Transpose[{steps,ConstantArray[0.0,nsteps+1]}];
colors=Table[ColorData[97][i],{i,nsteps+1}];
labels=Table[
Text[style[ToString[i]],curvepoints[[i]],{0,-1}],
{i,nsteps+1}
];
labels=Riffle[colors,labels];
axispoints=Riffle[colors,Point[{#}]&/@axispoints];
curvepoints=Riffle[colors,Point[{#}]&/@curvepoints];
rootpoints=Point[Transpose[{roots,{0,0}}]];
{
Plot[f4[x0][x],{x,xmin,xmax}
,PlotRange->{All,{-0.01,All}}
,PlotLegends->{f4[x0][x]}
,ImageSize->is
,Epilog->{
PointSize[Large],
curvepoints,
axispoints,
labels,
{Black,rootpoints}
}
],
steps
}
]//Column