odefun=@(t,y)[t-y];
tspan=[-1,5];
y0=1;
h=0.001;
Nh=round((tspan(2)-tspan(1))/h);
[tn,un]=beuler(odefun,tspan,y0,Nh);
figure(1)
plot(tn,un,'+')
xlabel('t')
ylabel('y(t)')

%%
% we define the exact solution
yex=@(t)[t-1+3*exp(-t-1)];
y=yex(tn);
err=max(abs(y-un))

hold on
plot(tn,y,'r')
legend('numerical sol','exact sol')






