%* routine to solve dy/dx=x+y using Euler steps

%* step size & number of intervals, N intervals, N+1 nodes
h=0.1;

%* independent variable
xi=0;
xf=5;
x=[xi:h:xf];			% this is the model domain
N=length(x)-1;			% number of steps, nodes=N+1

%* dependent variable
y0=1;					% initial condition
y=zeros(1,N+1);			% vector to store dependent variable
y(1)=y0;				% insert initial condition into y

%* numerical solution
for j=1:N
  y(j+1) = y(j) + h * (x(j) + y(j));
end

%* exact solution for x0=0, y0=1

g=2*exp(x) - x - 1;

%* compare results
figure(2)
clf
plot(x, y, 'b-')
hold on
plot(x, g, 'r-')
title('dy/dx = x + y')
legend('Euler single-step solution ', 'exact solution ')