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

/* User input comes first */
dydx: y/(2*x)-x/(2*y^2)   $/* dy/dx = f(x,y) */
x0: 3   $/* Initial x-value */
y0: 3   $/* 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
   (
   k1: h*ev( dydx, numer, [x=xn,y=yn] ),
   w: yn + k1,
   xn: xn + h,
   k2: h*ev( dydx, numer, [x=xn,y=w] ),
   yn: yn + ( k1 + k2 ) / 2,
   print("x", i, "=", xn, ", y", i, "=", yn)
   );