% variable time step RK23 scheme from homework
% with more than optimal number of feval function calls

function [s,n]=RK23hw(f,ab,y0,delta)
   yn=y0;
   h=(ab(2)-ab(1))/1000;
   t=ab(1);
   dab=delta/(ab(2)-ab(1));
   n=0;
   while t<ab(2)
       if t+h>ab(2)
          h=ab(2)-t;
       end
       n=n+1;
       xi1=yn;
       xi2=yn+h*1/2*feval(f,t,xi1);
       xi3=yn+h*(-1*feval(f,t,xi1)+2*feval(f,t+h/2,xi2));
       ynp1=yn+h*feval(f,t+h/2,xi2);
       xnp1=yn+h*((1/6)*feval(f,t,xi1) ...
               +(2/3)*feval(f,t+h/2,xi2) ...
               +(1/6)*feval(f,t+h,xi3));
       kappa=norm(ynp1-xnp1);
%       disp(sprintf('yn=%g t=%g kappa=%g h=%g h*dab=%g',...
%          yn,t,kappa,h,h*dab));
       if kappa<dab*h
         yn=ynp1;
         t=t+h;
         if kappa<0.1*dab*h
            h=h*2;
%            disp('increase h by doubling it');
         end
       else
         h=h/2;
%         disp('decrease h by halving it');
       end       
   end
   s=yn;
end