function approx = euler(h,y,t,f)

% Computes approximate solution to the given IVP for step size h 

% h is the step size
% y is the initial value
% t is the time interval: t = [t0 tf], where t0=initial time, tf=final time
% f is the right hand side of the system y' = f(t,y)

% Check input
if (h<=0)
    disp('Error: step size must be greater than 0.')
    return
end

t0 = t(1); tf = t(2);
n = (tf-t0)/h;                 % number of steps
approx = zeros(n+1,length(y)); % approximation of solution at each time step
approx(1,:) = y;               % initial value

% Implement forward Euler
for i=1:n,
    y = y + h*f(t0,y);   % approximation at step n, where time = t0 + n*h
    t0 = t0 + h;
    
    approx(i+1,:) = y;
end

