%4th order Taylor method for x'=1+x^2+t^3
function [ta,xa]=taylor4(t0,x0,tmax,npts)
n=npts-1;
h=(tmax-t0)/n;
ta(1)=t0;
xa(1)=x0;
for k=1:n
    t=ta(k);
    x=xa(k);
    xp=1+x^2+t^3;
    xpp=2*x*xp+3*t^2;
    xppp=2*x*xpp+2*xp*xp+6*t;
    xpppp=2*x*xppp+6*xp*xpp+6;
    x=x+h*(xp+0.5*h*(xpp+h*(xppp+0.25*h*xpppp)/3));
    t=t0+k*h;
    xa(k+1)=x;
    ta(k+1)=t;
end;