function [t, y] = myeuler(f, tinit, yinit, b, n)

% Euler approximation for initial value problem, 
% dy/dt = f(t, y), y(tinit) = yinit.
% Approximations y(1),..., y(n + 1) are calculated at
% the n + 1 points t(1),..., t(n + 1) in the interval
% [tinit, b]. The right-hand side of the differential
% equation is defined as an anonymous function f.

% Calculation of h from tinit, b, and n.

h = (b - tinit)/n;

% Initialize t and y as length n + 1 column vectors.

t = zeros(n + 1, 1);
y = zeros(n + 1, 1);

% Calculation of points t(i) and the corresponding
% approximate values y(i) from the Euler Method formula.

t(1) = tinit;
y(1) = yinit;
for i = 1:n
   t(i + 1) = t(i) + h;
   y(i + 1) = y(i) + h*f(t(i), y(i));
end