function [x, y] = myeuler(f, xinit, yinit, b, n)
   % Euler approximation for an IVP dy/dx = f(x,y)

     h = (b - xinit)/n;
     x = zeros(n+1,1);
     y = zeros(n+1,1);

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