clear all;
t0=0;
tN=2;
y0=1;
N=20;
h=(tN-t0)/N;
t=linspace(t0,tN,N+1);
exact=inline('t+exp(-t)', 't');
yexact=exact(t);
f=inline('-y+t+1', 't','y');
g=inline('y-t','t','y');
Ty(1)=y0;
for n=2:N+1
    Ty(n)=Ty(n-1) +h*f(t(n-1),Ty(n-1))+h^2/2*g(t(n-1),Ty(n-1));
end
plot(t,yexact,t,Ty,'+');
