function [tall, xall] = eulermethod_s19(f,tspan,x0,h)
% --------------------------------------------------------------
% Implements Euler's method for single and sytems of odes. 
% f -- an @ function (possibly vector-valued), 
% tspan -- a row vector with two elements, 
% x0 -- the initial condition vector,
% h -- the step size. 
%
% Returns a vector of time steps tall and approximations xall to 
% x(tall), where x(t) is the solution of x'=f(t,x) through (tspan(1),x0).
% ----------------------------------------------------------------

%-------------------------------------------------------------
% Define initial and final time t0 and tn, and number of steps
%-------------------------------------------------------------
t0=tspan(1);            
tn=tspan(2);
num_steps = (tn - t0)/h;

%-------------------------------------------------------------
% Set initial conditions tk,xk for Euler's method. 
% Initialize output vectors.  We set these up so that they
% have the same dimensions as the output from ode45: each
% row gives the value at a subsequent time value.
%-------------------------------------------------------------
tk = t0;
xk = x0;

% preallocate vectors for t and x
tall = ones( num_steps, 1 );
tall(1) = t0;

xall = ones( num_steps, length(x0) );
xall(1,:) = x0;

%-------------------------------------------------------------
% Loop for Euler's method. 
%-------------------------------------------------------------
for i=1:num_steps
    tall(i+1) = tall(i) + h;
    
    % recall: ' is the transpose operator in Matlab...
    xall(i+1,:) = xall(i,:) + h*f(tall(i), xall(i,:)).';
end

end

