function [x,y] = my_Euler2(fun,a,b,y0,N)
%  Euler method
%  Integration interval [a,b]
%  Initial condition y0
%  Number of (partition) points N+1

x = zeros(N+1,1) ;
y = zeros(N+1,1) ;

h = (b-a)/N ;
x(1) = a ;
y(1) = y0 ;
for n = 1:N
     x(n+1) = x(n) + h ;
     y(n+1) = y(n) + h*feval(fun,x(n),y(n)) ;
end    

end

