clear
f = @(t,y) t+y;
a = 0;
b = 1;

%

N = 5;
h = (b-a)/N;
t = linspace(a,b,N+1);
y(1) = 1;
for k = 2:N+1
    y(k) = y(k-1) + h*f(t(k-1),y(k-1));
end

%

plot(t,y,'or')

%

sol = @(x) 2*exp(x)- x - 1
hold on
fplot(sol,[0,1],'b')
hold off