function  [x, y] = rk4(f, ab, y0, N)
% Compute approximation of the solution of the initial value 
% problem  y' = f(x,y), y(a) = y0  on the interval ab = [a,b].
% Classical Runge-Kutta method with step length  h = (b-a)/N .

% Version 11.12.2003.  INCBOX

a = ab(1);  b = ab(2);  
h = (b - a)/N;  hh = h/2;  % full and half step length
x = linspace(a,b,N+1);     % grid
y = repmat(y0(:),1,N+1);   % initialize  y
for  n = 1 : N
  k1 = feval(f, x(n), y(:,n));
  k2 = feval(f, x(n)+hh, y(:,n)+hh*k1);
  k3 = feval(f, x(n)+hh, y(:,n)+hh*k2);
  k4 = feval(f, x(n+1), y(:,n)+h*k3);
  y(:,n+1) = y(:,n) + h/6*(k1 + 2*k2 + 2*k3 + k4);  
end
x = x';  y = y';           % return in standard format