function [u,t] = feuler( fun, f0, a, b, N )

    h=(b-a)/N;
    t = linspace(a,b,N+1);
    u(1) = f0;  
    
    for i=1:N,
        u(i+1) = u(i) + h * fun( t(i), u(i) );
    end

end






