%Assumes x0, y0, xfinal, and N have values before calling.
%The differential equation is defined below in the line that says yprime=.
%The result is 'xrecord' and 'yrecord', x and y values for the approximate
%solution.  Plot these with 'plot(xrecord,yrecord)'.

deltax=(xfinal-x0)/N;
xrecord=[x0];
yrecord=[y0];
x=x0;
y=y0;
for i=1:N
    yprime=x+y;
    x=x+deltax;
    y=y+yprime*deltax;
    xrecord=[xrecord x];
    yrecord=[yrecord y];
end