function [tout,yout] = myrk4(F,tspan,y0,h_arg,varargin)
%MYRK4  Classical fourth-order Runge-Kutta
%
%   Usage is the same as ODE23TX except the fourth argument
%   is a fixed step size h.
% 
%   MYRK4(F,TSPAN,Y0,H) with TSPAN = [T0 TFINAL] integrates the system
%   of differential equations y' = f(t,y) from t = T0 to t = TFINAL.  The
%   initial condition is y(T0) = Y0.  The input argument F is the
%   name of an M-file, or an inline function, or simply a character string,
%   defining f(t,y).  This function must have two input arguments, t and y,
%   and must return a column vector of the derivatives, y'.
%
%   With no output arguments, ODE23TX plots the emerging solution.
% 
%   With two output arguments, [T,Y] = MYRK4(...) returns a column 
%   vector T and an array Y where Y(:,k) is the solution at T(k).
%
%   More than four input arguments, MYRK4(F,TSPAN,Y0,RTOL,P1,P2,...),
%   are passed on to F, F(T,Y,P1,P2,...).
%
%   MYRK4(F,TSPAN,Y0,OPTS) where OPTS = ODESET('reltol',RTOL, ...
%   'abstol',ATOL,'outputfcn',@PLOTFUN) uses relative error RTOL instead
%   of 1.e-3, absolute error ATOL instead of 1.e-6, and calls PLOTFUN
%   instead of ODEPLOT after each successful step.
%
%   Example    
%      tspan = [0 2*pi];
%      y0 = [1 0]';
%      F = '[0 1; -1 0]*y';
%      myrk4(F,tspan,y0);
%
%   See also ODE23.
% 
% Written by:
% -- 
% John L. Weatherwax                2006-07-25
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

% Initialize variables.
%
h = 1.e-3; % a default step size
plotfun = @odeplot;
if nargin >= 4 & isnumeric(h_arg)
   h = h_arg;
end
t0 = tspan(1);
tfinal = tspan(2);
tdir = sign(tfinal - t0);
plotit = (nargout == 0);
t = t0;
y = y0(:);

% Make F callable by feval.

if ischar(F) & exist(F)~=2
   F = inline(F,'t','y');
elseif isa(F,'sym')
   F = inline(char(F),'t','y');
end 

% Initialize output.

if plotit
   feval(plotfun,tspan,y,'init');
else
   tout = t;
   yout = y.';
end

% The main loop
%
while t ~= tfinal
  
  % Evaluate the first slope:
  s1 = feval(F, t, y, varargin{:});

  % Stretch the step if t is close to tfinal:
  if 1.1*abs(h) >= abs(tfinal - t)
    h = tfinal - t;
  end
  
  % Evaluate the remaining slopes: 
  s2 = feval(F, t+h/2, y+(h/2)*s1, varargin{:});
  s3 = feval(F, t+h/2, y+(h/2)*s2, varargin{:});
  s4 = feval(F, t+h, y+h*s3, varargin{:});
   
  % Attempt a step.
  tnew = t + h;
  ynew = y + h*(s1 + 2*s2 + 2*s3 + s4)/6;
 
  tout(end+1,1) = tnew; 
  yout(end+1,:) = ynew.';
  
  t = tnew;
  y = ynew;
end

if plotit
   feval(plotfun,[],[],'done');
end