% RK4 (with more than optimal number of feval function calls)
function s=RK4(f,ab,y0,n)
   s=4;
   yn=y0;
   h=(ab(2)-ab(1))/n;
   t=ab(1);
   for k=1:n
       xi1=yn;
       xi2=yn+h*(1/2)*feval(f,t,xi1);
       xi3=yn+h*(1/2)*feval(f,t+h/2,xi2);
       xi4=yn+h*feval(f,t+h/2,xi3);
       yn=yn+h*((1/6)*feval(f,t,xi1) ...
               +(1/3)*feval(f,t+h/2,xi2) ...
               +(1/3)*feval(f,t+h/2,xi3) ...
               +(1/6)*feval(f,t+h,xi4));
       t=ab(1)+k*h;
   end
   s=yn;
end