function [t,y] = rk4_sys(f, lh, ic, stepsize)
%Usage: [t,y] = rk4_sys(f, lh, ic, stepsize)
%RK4 Method implementation for ODE system y' = f(t,y) 
%hence y can be a vector. Here lh is the interal lh[1] <= t <= lh[2] over 
%which the solution is computed and y(lh[1]) = ic, where ic is the vector
%of initial conditions at time t = lh[1]. The step size is in argument 
%"stepsize". If the number of steps from initial to final time is not an 
%integer, the last step is incremental of proper length. Solution times 
%are in t, solution values in y.

%Prep for marching out in time
finalT = lh(2);
N = floor((finalT-lh(1))/stepsize); %Round down
T1 = lh(1)+N*stepsize; %Corresponding time
Nf = N+1; %Number of steps we'll take, unless T1 = finalT
if abs(T1-finalT) < 1.0e-8 %To some tolerance
    Nf = N;
end

%Make that ic is a row vector.
ic = ic(:)';

%Set up arrays for solution return
N = length(ic);
t = zeros(Nf+1,1);
y = zeros(Nf+1,N);
t(1) = lh(1); %Initial time
y(1,:) = ic; %Initial value

%March out in time
for i=2:Nf
    t(i) = t(i-1) + stepsize; %Update time
    
    %Extrapolate solution using RK4 formulae
    told = t(i-1); yold = y(i-1,:);
    m1 = f(told,yold); m1 = m1(:)';
    m2 = f(told+0.5*stepsize,yold+0.5*stepsize*m1); m2 = m2(:)';
    m3 = f(told+0.5*stepsize,yold+0.5*stepsize*m2); m3 = m3(:)';
    m4 = f(told+stepsize,yold+stepsize*m3); m4 = m4(:)';
    m = (m1+2*m2+2*m3+m4)/6.0;
    y(i,:) = yold + stepsize*m;
end

%Last step, might be special/incremental.
t(Nf+1) = finalT; %Update time
stepsize = finalT - t(Nf);
told = t(Nf); yold = y(Nf,:);
m1 = f(told,yold); m1 = m1(:)';
m2 = f(told+0.5*stepsize,yold+0.5*stepsize*m1); m2 = m2(:)';
m3 = f(told+0.5*stepsize,yold+0.5*stepsize*m2); m3 = m3(:)';
m4 = f(told+stepsize,yold+stepsize*m3); m4 = m4(:)';
m = (m1+2*m2+2*m3+m4)/6.0;
y(Nf+1,:) = yold + stepsize*m;