function [tout, yout] = ode_Runge(FunFcn, t, y0, steps)
%ODE_RUNGE
%       ode_runge(FunFcn, t, y0, steps)
%
%	Integrate a system of ordinary differential equations using
%	4th order Runge-Kutta formula.
%
% INPUT:
% F     - String containing name of user-supplied problem description.
%         Call: yprime = fun(t,y) where F = 'fun'.
%         T      - Vector of times (scalar).
%         y      - Solution column-vector.
%         yprime - Returned derivative column-vector; yprime(i) = dy(i)/dt.
% t(1)  - Initial value of t.
% y0    - Initial value column-vector.
% steps - steps to take between given output times
%
% OUTPUT:
% T  - Returned integration time points (column-vector).
% Y  - Returned solution, one solution column-vector per tout-value.
%
% The result can be displayed by: plot(tout, yout).

% Initialization
h = (t(2)-t(1)) /steps;
y = y0(:);
[n,m]=size(y);
if (n<m) y = y'; end
yout = y.';
tout= t(:);

% The main loop
for i = 2:length(t)
   h = (t(i)-t(i-1))/steps;
   tau = t(i-1);
   for j=1:steps
      % Compute the slopes
      s1 = feval(FunFcn, tau, y); s1 = s1(:);
      s2 = feval(FunFcn, tau+h/2, y+h*s1/2); s2 = s2(:);
      s3 = feval(FunFcn, tau+h/2, y+h*s2/2); s3 = s3(:);
      s4 = feval(FunFcn, tau+h, y+h*s3); s4 = s4(:);
      tau = tau + h;
      y = y + h*(s1 + 2*s2+ 2*s3 + s4)/6;
    end
    yout = [yout; y.'];
end
