Newton's Method for approximating a solution of f(x)=0.

/* User input comes first */
f: cos(x)-x   $/* Function f(x) */
x0: 1   $/* Initial estimate */
N: 5   $/* Number of iterations */
/* End of user input */

fp: diff(f,x,1)$
xn: float(x0)$
print("x", 0, "=", xn)$
for i: 1 thru N do
   (
   xn: xn - at(f,x=xn) / at(fp,x=xn),
   print("x", i, "=", xn)
   );

Riemann Sums to approximate the definite integral of f(x) on [a,b]

/* User input comes first */
f(x):= x^2+2*x   $/* Define function */
a: 0   $/* Left endpoint */
b: 1   $/* Right endpoint */
N: 20   $/* Number of subintervals */
opt: 1   $/* 0 for left endpts, 1 for midpts, 2 for right endpts */
/* End of user input */

h: (b-a)/N$
ans1: 'sum( f(a+(2*i+opt)*h/2)*h, i, 0, N-1 )$
ans2: ev(%, sum)$
ans3: float(%)$
print( ans1, " = ", ans2, " = ", ans3 )$

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*at(dydx,[x=xn,y=yn]),
   xn: xn + h,
   print("x", i, "=", xn, ", y", i, "=", yn)
   );

4th-order Runge-Kutta Method for numerically approximating solutions of the IVP:
dy/dx=f(x,y), y(x0)=y0

/* User input comes first */
dydx: -1.2*y-9.8   $/* dy/dx = f(x,y) */
x0: 0   $/* Initial x-value */
y0: 20   $/* Initial y-value */
Xend: 2   $/* Ending x-value */
steps: 20   $/* Number of steps --- stepsize will be (Xend - x0) / steps */
/* End of user input */

xx: float(x0)$ yy: float(y0)$
h: ( float(Xend) - float(x0) ) / steps$
print("Step size = ", h)$
print("x",0,"=", xx, ", y",0,"=", yy)$
for i: 1 thru steps do
   (
   k1: h*at(dydx,[x=xx,y=yy]),
   ww: yy + k1/2,
   k2: h * at(dydx,[x=xx+h/2,y=ww]),
   ww: yy + k2/2,
   k3: h * at(dydx,[x=xx+h/2,y=ww]),
   ww: yy + k3,
   k4: h * at(dydx,[x=xx+h,y=ww]),
   yy: yy + ( k1 + 2*k2 + 2*k3 + k4 ) / 6.0,
   xx: xx + h,
   print("x", i, "=", xx, ", y", i, "=", yy)
   );