ODEs
(*Euler's method *)

ODE: y’ = f(x,y).
Approximate: (ynew-yold)/h = f(xold, yold)
domain: (x0,xn)

numpts= 50; (* how many evals *)
x0=0;
xn=2;
h=(xn-x0)/numpts//N;
f[x_,y_]:=y;
(*initialze array: *)
eularray = Table[1,{i,1,numpts}];	
   
Do[ (* fill eularray *)
(
xold= (i-1)*h;
xnew=i*h;
yold=eularray[[i-1]];
ynew=yold + (f[xold,yold])*h;
eularray[[i]]=ynew;
),
{i,2,numpts}
]

exactplot=Plot[E^x,{x,x0,xn},PlotStyle->{GrayLevel[0.4],Thickness[0.01]}] ; 
(* example E^x *)

numpts1= 100;
h1=(xn-x0)/numpts1//N;
f[x_,y_]:=y;

Print["last Euler values: ",eularray[[numpts]],"  ",eularray1[[numpts1]]];
Print["exact last value: ", E^xn//N];




Show[{ exactplot}]

ListPlot[{eularray, eularray1},
PlotStyle->Directive[PointSize[Medium]]]
(* note different scales in x *)