function [h,t,u] = expeuler(fun, t0, y0, N, a)

h=a/N;
t=t0+(0:N)*h;
u=zeros(length(y0),N+1); % List initialization
u(:,1)=y0;
for i=1:N % Start the Explicit Euler procedure
    ui= u(:,i);
    ti=t(i);
    u(:,i+1)=ui+h*fun(ti, ui); % y is being overwritten by the new iterate
end
end