function [tout,yout,dt,success]=step(time,y,dt,tol,tmax,p)
global diff;
pm=p-1;
[t1,y1]=rk4(time,y,dt/2);
[t2,y2]=rk4(t1,y1,dt/2);
[t3,y3]=rk4(time,y,dt);
%diff=sqrt((y3(1)-y2(1))^2+(y3(2)-y2(2))^2);
diff=y3(1)-y2(1);
errorestimate=abs(diff)/(1-2^(-pm));
dtopt=((tol/tmax)*(dt^p)/errorestimate)^(1/pm);
if dtopt > dt
  tout=t3;
  yout=y2;
%  yout=2*y2-y3
  yout=16*y2/15-y3/15;
%  yout=y2;
  success=1;
  dt=min(.9*dtopt,2*dt);
else
  tout=time;
  yout=y;
  [time,dt,dtopt]
  dt=.9*dtopt;
  success=0;
end

 