% Basic Runge-Kutta 4th order ODE solver 
function [yout, xout] = RK4(dy, y0, xspan)
  OneDim=0;
  if min(size(y0))==1
    OneDim=1;
    y = y0(:);
  else
    y = y0;
  end
  if length(xspan) == 2
    xspan = linspace(xspan(1),xspan(2),101);
  end

  xout = xspan;
  yout = y;

  n = 1;
  N = length(xspan);
  while (n < N)
    x=xspan(n);
    dx = diff(xspan(n:n+1));

    % Compute the slopes
    k1 = dy(y,x);
    k2 = dy(y + dx*k1/2,x + dx/2);
    k3 = dy(y + dx*k2/2,x + dx/2);
    k4 = dy(y + dx*k3,x + dx);
    y = y + dx*(k1 + 2*(k2 + k3) +k4)/6;

    if OneDim
      yout(:,end+1) = y;
    else
      yout(:,:,end+1) = y;
    end
    n=n+1;
  end
end
