function [y,dydx,istop] = rungeKutta_JMS( deriv, y0, x0, dx, xmax, stopCheck )

% Solves first-order differential equation dy/dx=deriv(x,y), using
% fourth-order Runge-Kutta method, where y is a column vector and 
% x is a scalar variable.
% J.M.Soler, May 2019
% 
% Input:
%   deriv     : function returning the derivative as dy/dx=deriv(x,y)
%   y0(ny,1)  : initial y values
%   x0        : initial x value
%   dx        : integration step
%   xmax      : final x value
%   stopCheck : stopping function called as istop=stopCheck(x,y)
%               where istop~=0 will stop the trajectory (optional)
% Output:
%   y(ny,nx)    : y values for x=(x0:dx:xmax)
%   dydx(ny,nx) : derivative dy/dx at x=(x0:dx:xmax)
%   istop       : code returned by stopCheck, or zero, if not present

% Algorithm:
% The fourth-order Runge-Kutta iteration is
%   z1 = y'(x,y),           where y=y(x)  and  y'=dy/dx,
%   z2 = y'(x+dx/2,y+z1*dx/2),
%   z3 = y'(x+dx/2,y+z2*dx/2),
%   z4 = y'(x+dx,y+z3*dx).
%   y(x+dx) = y(x) + z1*dx/6 + z2*dx/3 + z3*dx/3 + z4*dx/6
% Ref: W.H.Press et al, Numerical Recipes, Cambridge U.P.

% Find array sizes
nx = numel(x0:dx:xmax);                 % number of x values
ny = numel(y0);                         % size of vector y

% Integrate trajectory
istop = 0;                              % default value
x = x0;                                 % initial values
y(:,1) = reshape(y0,ny,1);
dydx(:,1) = deriv(x0,y(:,1));
for ix = 2:nx
  y1 = y(:,ix-1);
  z1 = dydx(:,ix-1);
  z2 = deriv(x+dx/2,y1+z1*dx/2);
  z3 = deriv(x+dx/2,y1+z2*dx/2);
  z4 = deriv(x+dx,y1+z3*dx);
  x = x+dx;
  y(:,ix) = y1 + z1*dx/6 + z2*dx/3 + z3*dx/3 + z4*dx/6;
  dydx(:,ix) = deriv(x,y(:,ix));
  if nargin>5                           % check stop condition
    istop = stopCheck(x,y(:,ix));
    if istop~=0                         % stop trajectory
      break                             % exit ix loop
    end
  end
end

end % function rungeKutta
