euler_method(f, y0, a, b, h):= block(
  [t: a, y: y0, tg: [a], yg: [y0]],
  unless t>=b do (
    t: t + h,
    y: y + f(t, y)*h,
    tg: endcons(t, tg),
    yg: endcons(y, yg)
    ),
  [tg, yg]
  );

/* initial temperature */
T0: 100;

/* environment of temperature */
Tr: 20;

/* the cooling constant */
k: 0.07;

/* end of integration */
tmax: 100;

/* analytical solution */
Tref(t):= Tr + (T0 - Tr)*exp(-k*t);

/* cooling rate */
dT(t, T):= -k*(T-Tr);

/* get numerical solution */
h: 10;
[tg, yg]: euler_method('dT, T0, 0, tmax,  h);

/* plot analytical and numerical solution */
plot2d([Tref, [discrete, tg, yg]], ['t, 0, tmax],
  [legend, "analytical", concat("h = ", h)],
  [xlabel, "t / seconds"],
  [ylabel, "Temperature / C"]);