%
%   EULER.M   Use Euler Method for Numerical Solution of IVPs
%
%   This function file is modeled after the example in by Edwards and 
%   Penny (see pg. 105).  However it is designed to handle vector functions.
%
%   The inputs to the function are:
%      fxy = string variable with the name of the function file containing f(x,y)
%      xo,xf = initial and final values of the independent variable (scalars)
%      yo = initial value of dependent variable at xo (column vector)
%      N = number of intervals to use between xo and xf
%
%   The outputs to the function are:
%      X = vector containing values of the independent variable
%      Y = the estimated dependent variable at each value of the independent variable
%       --> this variable is a vector if only one equation is solved
%       --> it is a matrix [y1(x) y2(x) ... ] for multiple equations 
%
%   File prepared by J. R. White, UMass-Lowell (March 1998)
%

      function [X,Y] = euler(fxy,xo,xf,yo,N)
%
%   compute step size and size of output variables
      if N < 2   N = 2;   end  % set minimum number for N
      h = (xf-xo)/N;           % step size
      X = zeros(N+1,1);        % initialize independent variable
      M = max(size(yo));       % number of equations (number of columns of Y matrix)
      Y = zeros(N+1,M);        % initialize dependent variables
%
%   set initial conditions
      x = xo;   X(1) = x;   y = yo;   Y(1,:) = y';
%
%   begin computational loop
      for i = 1:N
        k1 = h*feval(fxy,x,y);   % evaluate function
        y = y + k1;              % increment dependent variable
        x = x + h;               % increment independent variable
        X(i+1) = x;              % store current x in X vector
        Y(i+1,:) = y';           % store current y in Y matrix
      end
      disp([X Y])                % This line added by Steve Pennell, Oct. 2010  
%
%   end of function
%