>function dgleuler (f,x,y0) ...
$  y=zeros(size(x)); y[1]=y0;
$  for i=2 to cols(y);
$  y[i]=y[i-1]+f(x[i-1],y[i-1])*(x[i]-x[i-1]);
$  end;
$  return y;
$endfunction
>function f(x,y) := -k*(y-TR)
>k=0.07; TR=20; TS=100;
>x=0:1:100; dgleuler("f",x,TS)[-1]
 20.0564137335
>x=0:2:100; dgleuler("f",x,TS)[-1]
 20.0424631834
>TR+(TS-TR)*exp(-k*TS)
 20.0729505572
>x=0:5:100; plot2d(x,dgleuler("f",x,TS)); ...
>  plot2d(x,TR+(TS-TR)*exp(-k*x),>add,color=red);
>ode("f",x,TS)[-1] // Euler default solver LSODA
 20.0729505568
>adaptiverunge("f",x,TS)[-1] // Adaptive Runge Method
 20.0729505572