Taylor Methods applied to the motion of a projectile,
p. 344, Computer problem #1

Taylor-4 with step-size control

Clear[f1,f2,f3,f4,g,r,v,v1,v2,v3,x]
(*
This is the function, representing the rate of change of the velocity of a projectile fired from the surface of the earth with respect to height, where x is the height above the earth. 

As an aside, this example I worked on while listening to Nobel prize-winning physicist Richard Feynman talk about the Law of Universal Gravitation, https://www.youtube.com/watch?v=JzWzLyGuPRY
 *)
f1=-g r^2/(v[x] (r+x)^2)
(* now we compute some partial derivatives, implicitly: *)
f2=D[f1,x];
f3=D[f2,x];
f4=D[f3,x];
(* and now we change the formal derivatives into atomic symbols *)
f1=f1/.v[x]->v
f2=f2/.v[x]->v;
f2=f2/.v'[x]->v1
f3=f3/.v[x]->v;
f3=f3/.v'[x]->v1;
f3=f3/.v''[x]->v2
f4=f4/.v[x]->v;
f4=f4/.v'[x]->v1;
f4=f4/.v''[x]->v2;
f4=f4/.v'''[x]->v3

Some physical parameters: radius of the Earth r, and the acceleration due to gravity g

r=6378.388; (* in km *)
g=0.00980621;  (* in km/s^2 *)

Here’s Taylor’s 4 method:

taylor4[vi_,xf_]:=
Module[
{xs,vs,ltol=0.0001,maxsteps=100,k=0,h,xnew},
xs={x=0.};
vs={v=vi};
While[x<xf&&v>0.&&k<maxsteps, (* height is less than xf, the velocity is positive -- we're not falling back, and we've not taken more than the allowed number of step: *)
(* we compute the higher derivatives needed for Taylor-4: *)
v1=N[f1];
v2=N[f2];
v3=N[f3];
v4=N[f4];
(* now we introduce some automatic step-size control, based on the error: *)
h=(4! ltol/Abs[v4])^0.25;
(* check that we don't step past the final value of x: *)
xnew=Min[x+h,xf];
(* adjust h if it happens to be too big, stepping beyond the final x value *)
h=xnew-x;
(* Here we evaluate our next estimate for the position of our projectile: note the use of Horner's method to evaluate the polynomial of Taylor-4: 
*)
v=v+h (v1+h/2. (v2+h/3. (v3+h/4. v4)));
x=xnew;
AppendTo[xs,x];
AppendTo[vs,v];
k=k+1];
{xs,vs}]

Now we’ll look at the two examples suggested in the text:

(* here we're starting with the escape velocity *)
vi=Sqrt[2. g r] (* in km/s *)
xf=10.^6;
t4=taylor4[vi,xf] ;
(* This equatio ncan be solved by separation of variables: *)
vTheoretical[x_]:=Sqrt[2g*r^2/(r+x) +vi^2-2g*r]
p0=Plot[vTheoretical[x],{x,0,xf},PlotRange->{{0,xf},{0,vi}}];
p1=ListPlot[Transpose[t4],PlotMarkers->{Automatic,Tiny}];
Show[p0,p1]
(*The foregoing computation is a striking illustration of the power of automatic stepsize control.*)

(* here we're starting with an initial velocity just below escape velocity: *)
vi=11. (* in km/s *)
p0=Plot[vTheoretical[x],{x,0,200000},PlotRange->{{0,200000},{0,vi}}]
t4=taylor4[vi,xf];
p1=ListPlot[Transpose[t4],PlotMarkers->{Automatic,Medium}];
Show[p0,p1]