function [X,Y] = euler(f,span,y,n)
h = (span(2)-span(1))/n;  % step size
x = span(1);              % initial x
X = x;                    % initial x
Y = y;                    % initial y
for i = 1:n               % begin loop
    y = y + h*f(x,y);     % Euler iteration
    x = x + h;            % new x
    X = [X;x];            % update x-column
    Y = [Y;y];            % update y-column
end