t=0; y=1; dt=0.01; T=[t,y];
while t<0.8
  y1=t^2+y^2;
  y2=2*t+2*y*y1;
  y3=2+2*y1^2+2*y*y2;
  y4=6*y1*y2+2*y*y3;
  y=y+y1*dt+y2/2*dt^2+y3/6*dt^3+y4/24*dt^4;
  t=t+dt;
  T=[T;t,y];
end;
[t,Y]=ode45(@(t,y)t^2+y^2,[0:dt:0.8],1);
norm(T(:,2)-Y)

