function [tout, yout] = rk4(FunFcn, tspan, y0, ssize)
%[TOUT,YOUT] = rk4(ODEFUN,TSPAN,Y0,OPTIONS)
% Solve non-stiff differential equations, medium order method.
% with TSPAN = [T0 TFINAL] integrates 
% the system of differential equations y' = f(t,y) from time T0 to TFINAL 
% with initial conditions Y0. ODEFUN is a function handle. For a scalar T
% and a vector Y, ODEFUN(T,Y) must return a column vector corresponding 
% to f(t,y). Each row in the solution array YOUT corresponds to a time 
% returned in the column vector TOUT.  To obtain solutions at specific 
% times T0,T1,...,TFINAL (all increasing or all decreasing), use TSPAN = 
% [T0 T1 ... TFINAL].  
%exemple:
%[t,y]=rk4(@vdp1,[0 20],[2 0]);   
%plot(t,y(:,1));

% Initialization
t0 = tspan(1);
tfinal = tspan(2);
pm = sign(tfinal - t0); % Which way are we computing?
if nargin < 4, ssize = (tfinal - t0)/100; end
if ssize < 0, ssize = -ssize; end
h = pm*ssize;
t = t0;
y = y0(:);
% We need to compute the number of steps.
dt = abs(tfinal - t0);
N = floor(dt/ssize) + 1;
if (N-1)*ssize < dt
N = N + 1;
end
% Initialize the output.
tout = zeros(N,1);
tout(1) = t;
yout = zeros(N,size(y,1));
yout(1,:) = y.';
k = 1;
% The main loop
while (k < N)
if pm*(t + h - tfinal) > 0
h = tfinal - t;
tout(k+1) = tfinal;
else
tout(k+1) = t0 +k*h;
end
k = k + 1;
% Compute the slopes
s1 = feval(FunFcn, t, y); s1 = s1(:);
s2 = feval(FunFcn, t + h/2, y + h*s1/2); s2 = s2(:);
s3 = feval(FunFcn, t + h/2, y + h*s2/2); s3 = s3(:);
s4 = feval(FunFcn, t + h, y + h*s3); s4 = s4(:);
y = y + h*(s1 + 2*s2 + 2*s3 +s4)/6;
t = tout(k);
yout(k,:) = y.';
end;