%This is eul.m

function [tout, yout] = eul(FunFcn, t0, tfinal, y0, ssize)

% usage:  [tout, yout] = eul(FunFcn, t0, tfinal, y0, ssize)
%
% EUL	Integrates a system of ordinary differential equations using
%	Euler's method.  See also ODE45 and
%	ODEDEMO.M.
%	[T,Y] = EUL('yprime', T0, Tfinal, Y0) integrates the system
%	of ordinary differential equations described by the M-file
%	YPRIME.M over the interval T0 to Tf and using initial
%	conditions Y0.
%	[T, Y] = EUL(F, T0, Tfinal, Y0, SSIZE) uses step size SSIZE
%
% INPUT:
% F     - String containing name of user-supplied problem description.
%         Call: yprime = fun(t,y) where F = 'fun'.
%         t      - Time (scalar).
%         y      - Solution vector.
%         yprime - Returned derivative vector; yprime(i) = dy(i)/dt.
% t0    - Initial value of t.
% tfinal- Final value of t.
% y0    - Initial value vector.
% ssize - The step size to be used. (Default: ssize = (tfinal - t0)/100).
%
% OUTPUT:
% T  - Returned integration time points (column-vector).
% Y  - Returned solution, one solution row-vector per tout-value.
%
% The result can be displayed by: plot(tout, yout).

% Written by John C. Polking, Rice University
% Last modifed February 3, 1995


% Initialization

if (nargin < 5), ssize = (tfinal - t0)/100; end

h = ssize;
t = t0;
y = y0(:);
tout = t;
yout = y.';

% The main loop
while (t < tfinal)
      if t + h > tfinal, h = tfinal - t; end

      % Compute the slope
      s1 = feval(FunFcn, t, y); s1 = s1(:); % s1=f(t(k),y(k))

      t = t + h;
      y = y + h*s1;		% y(k+1) = y(k) + h*f(t(k),y(k))
      tout = [tout; t];
      yout = [yout; y.'];

end;