Off[General::spell1];
Clear[TaylorM ,f0,a,b,x0,n];
TaylorM[f0_,a_,b_,x0_,n_] :=
Module[{d1,d2,d3,d4,f,h,j,T, X,X2,X3,X4},
  Set@@{f[t_,x_],f0};
  X'[t] = f[t,X[t]];
  X2=D[X[t],{t,2}]; 
  X3=D[X[t],{t,3}]; 
  X4=D[X[t],{t,4}]; 
  Set@@{d1[t_,x_],f[t,x]};
  Set@@{d2[t_,x_], X2/. X[t]->x };
  Set@@{d3[t_,x_], X3/. X[t]->x };
  Set@@{d4[t_,x_], X4/. X[t]->x };
  h = (b - a)/n //N;
  T = Table[0,{n+1}];
  X = Table[0,{n+1}];
  T[[1]] = a //N;
  X[[1]] = x0 //N;
  For[j=1, j<=n,j++,
    tj = T[[j]] //N;
    xj = X[[j]] //N;
    T[[j+1]] = a  + h j //N;
    X[[j+1]] = xj + h(d1[tj,xj] + 
     h(d2[tj,xj]/2 + h(d3[tj,xj]/6 +
     h d4[tj,xj]/24))); ];
  Return[Transpose[{T,X}]] ];
On[General::spell1];

Clear[f,t,x];   
f[t_,x_] = Cos[t] - Sin[x] +t^2;

TaylorM[f[t,x],-1,1,3,10]

Clear[pts];
pts = TaylorM[f[t,x],-1,1,3,100];
ListPlot[pts,PlotJoined->True, AspectRatio -> 0.5];

Clear[ss];
ss=NDSolve[{x'[t]==f[t,x[t]],x[-1]==3},x[t],{t,-1,1}];
Plot[Evaluate[x[t]/.ss],{t,-1,1},AspectRatio -> 0.5];

Clear[xx];
xx[t_]= x[t] /. First[ss];
Error = Table[pts[[i,2]]-xx[-1+ (i-1) /50], {i,1,10}]

Clear[f,t,x];
f[t_,x_] = (t *x - x^2) / t^2;

Clear[pts];
pts = TaylorM[f[t,x],1,3,2,100 ];
ListPlot[pts,PlotJoined->True, AspectRatio -> 1];

Clear[ss,xx];
ss=DSolve[{x'[t]==f[t,x[t]],x[1]==2},
   x[t],t];
xx[t_]= x[t] /. First[ss]

Plot[xx[t], {t, 1,3},AspectRatio -> 1]

Clear[xx];
xx[t_]= x[t] /. First[ss];
Error = Table[pts[[i,2]]-xx[1+ (i-1) /50], {i,1,10}]//N