Clear["Global`*"];

f[t_,x_]:=Sin[x] (t^3)+Tanh[t];                                      (* Define your f[t,x] here,  where the differential equation is  x'=f[t,x]   *)

f2[t_,x_]:=D[f[t,x],t]+D[f[t,x],x] f[t,x];            (*  This is the first derivative of f  (or the 2nd derivative of x)  *)

f3[t_,x_]:=D[f2[t,x],t]+D[f2[t,x],x] f[t,x];          (*  This is the 2nd derivative of f  (or the 3rd derivative of x)  *)

f4[t_,x_]:=D[f3[t,x],t]+D[f3[t,x],x] f[t,x];           (*  This is the 3rd derivative of f  (or the 4th derivative of x)  *)

xtph=x0+h f[t0,x0]+(1/2) (h^2) f2[t0,x0]+(1/(3!)) (h^3)f3[t0,x0]+(1/(4!)) (h^4)f4[t0,x0]                   (*  So the approximation to x at  t0+h  is  xtph  *)

xtph=Collect[Simplify[xtph],h]                                               (*  Here we simplify the expression above  a bit  *)

xtph/.{t0->0,x0->1,h->1/2}                                                
(*  Now, to obtain our numerical approximation to x at t=t0+h, 
 we substitute t0->0, x0->1,    and h->1/2  into the expression xtph *)

taylorSeriesApprox=N[%,10]                                                             (*  Get a numerical value for the above approximation  *)

tfinal = 3;   (*  we are interested in a solution valid over {0,3.0}  *)

(*  so, for comparison, we can get a very accurate solution using the 
built-in Mathematica function NDSolve the over the range {0, tfinal}  *)

solution=NDSolve[{x'[t]==Sin[x[t]] (t^3)+Tanh[t],x[0]==1},x,{t,0,tfinal}]                  (*  Now, for comparison, we could also use the built-in Mathematica function NDSolve *)
(*  to obtain a result which is good to machine precision  *)

(* and we can now plot this "true solution" returned from NDSolve  *)

m1plot=Plot[x[t]/.solution,{t,0,tfinal},PlotStyle->Hue[0]]

xTrue=x[t]/.solution/.t->1/2                                    (*  We evaluate our "true solution"  at  t=1/2  so we can compare with the Taylor Series approximation  *)

error=taylorSeriesApprox-xTrue            (*  and we find that the error is quite small  *)

(*  and we now plot our Taylor Series approximation  *)
m2plot=Plot[xtph/.{t0->0,x0->1,h->t},{t,0,tfinal},PlotStyle->Hue[0.7]]

Show[m1plot,m2plot,PlotRange->{0,10}]                                                                                         (*  and we compare the TS with the "true" solution   *)

(*  Notice how the Taylor Series is very good near the initial point (t0)  and deteriorates as we get farther from this point.  So, we should take a number of small steps in order to solve the differential equation accurately over the range of the plot shown above   *)

h=1;
t0=0.0;
x0=1.0;
list={};
Print[0,"   ",t0,"   ",x0];
(*  initialize a list for later plotting  *)
list=Append[list,{t0,x0}];
(*  initialize a list for later plotting  *)
(*  now, let's take tfinal/h steps  *)
For[i=0,i<tfinal/h,i++,
x0=xtph;
t0=t0+h;
Print[i,"   ",t0,"   ",x0];
(*  Now, add the new point to the list for plotting  *)
list=Append[list,{t0,x0}];   
];

lp1=ListPlot[list,PlotStyle->Darker[Green],PlotRange->{{-0.1,3.1},{-12.3,4}}]

h=0.5;
t0=0.0;
x0=1.0;
list={};
Print[0,"   ",t0,"   ",x0];
(*  initialize a list for later plotting  *)
list=Append[list,{t0,x0}];
Print[0,"   ",t0,"   ",x0];
(*  initialize a list for later plotting  *)
(*  now, let's take tfinal/h steps  *)
For[i=0,i<tfinal/h,i++,
x0=xtph;
t0=t0+h;
Print[i,"   ",t0,"   ",x0];
(*  Now, add the new point to the list for plotting  *)
list=Append[list,{t0,x0}];   
];

(*  we can just plot the values of {t0,x0} at each step   *)

lp2=ListPlot[list,PlotStyle->Blue,PlotRange->{{-0.1,3.1},{-12.2,4}}]

(*  Now, we can compare our two numerical solutions with the "true" solution   *)
(*  The results for the h=1 case are in green  *)
(*  The results for the h=0.5 case are in blue  *)
(*  The "true" solution is in red  *)

Show[m1plot,lp2,lp1,PlotRange->{{-0.1,3.1},{-12.2,4}}]

(*  and we can now observe that by taking small steps, we can follow the "true" solution much better  *)