function [x,y] = EulerFixH(dfun,xspan,y0,nsteps)
% [x,y] = EulerFixH(dfun,xspan,y0,nsteps)
%
% Integrates ODE system from x0 to xfinal by nsteps Euler steps
% Inputs:
%   dfun    String that names a function M-file for derivatives
%           in y'=f(x,y) so that feval(dfun,x,y) evaluates f(x,y)
%   xspan   2-component vector: interval on which to compute sol'n
%   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 (Euler approximation) on the mesh x
%
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
% 
% Euler's method loop - march across the interval.
% 
for n = 1:nsteps                   
  y(n+1,:) = [ y(n,:)' + h*feval(dfun,x(n),y(n,:)') ]';
end