function [x,y] = RK2FixH(dfun,xspan,y0,nsteps)
% [x,y] = RK2FixH(dfun,xspan,y0,nsteps)
% Integrates ODE on interval xspan by nsteps
% steps of a second-order Runge-Kutta method.
% Inputs:
%   dfun    string naming a function M-file for derivatives
%   xspan   integration interval
%   y0      initial value(s) of dependent variable(s)
%   nsteps  number of steps to take
% Outputs:
%   x       mesh of nsteps+1 equally spaced pts from x0 to xfinal
%   y       solution of DE (RK2 approximation) on the mesh x
%
c2 = 2/3;                            % Determines Runge-Kutta 
b2 = 1/(2*c2); b1 = 1-b2;            %  coefficients.
y0 = y0(:);                          % Column vector
x = linspace(xspan(1),xspan(2),nsteps+1)';  % Create mesh
[neq one] =size(y0);                 % Determine # of dep. vars.
y = [y0'; zeros(nsteps,neq)];        % Allocate space for solution
h = diff(xspan)/nsteps;              % Compute stepsize
for n = 1:nsteps                     % Step across the interval
  k1 = h*feval(dfun,x(n),y(n,:)');
  k2 = h*feval(dfun,x(n)+c2*h,...
               y(n,:)'+c2*k1);
  y(n+1,:) = [ y(n,:)' + b1*k1 + b2*k2 ]';
end