Euler's Method for numerically approximating solutions of the IVP:
dy/dx=f(x,y), y(x0)=y0

/* User input comes first */
dydx: x*y   $/* dy/dx = f(x,y) */
x0: 0   $/* Initial x-value */
y0: 1   $/* Initial y-value */
h: 0.1   $/* Step size */
steps: 10   $/* Number of steps -- Ending x-value will be x0 + steps*h */
/* End of user input */

xn: float(x0)$ yn: float(y0)$
print("x", 0, "=", xn, ", y", 0, "=", yn)$
for i: 1 thru steps do
   (
   yn: yn + h*ev( dydx,numer, [x=xn,y=yn] ),
   xn: xn + h,
   print("x", i, "=", xn, ", y", i, "=", yn)
   );