# Solving differential equations

# See Handbook of Mathematics and Computational Science,
# John W.Harris, Horst Stocker
SetHelp ("EulersMethod", "equation_solving",
         "Use classical Euler's method to numerically solve y'=f(x,y) for initial x0,y0 going to x1 with n increments, returns y at x1")
function EulersMethod(f,x0,y0,x1,n) = (
	# Note we can't check the 2 arguments, FIXME
	if not IsFunction(f) then
		(error("EulersMethod: f must be a function of two arguments");bailout)
	else if not IsValue(x0) or not IsValue(y0) or not IsValue(x1) then
		(error("EulersMethod: x0, y0 and x1 must be numbers");bailout)
	else if not IsPositiveInteger(n) then
		(error("EulersMethod: n must be a positive integer");bailout);
	h := float(x1-x0) / n;
	x := x0;
	y := y0;
	for k = 0 to n do (
		y := y + h*f(x,y);
		x := x + h
	);
	y
)
protect ("EulersMethod")

# See Handbook of Mathematics and Computational Science,
# John W.Harris, Horst Stocker
SetHelp ("RungeKutta", "equation_solving",
         "Use classical non-adaptive Runge-Kutta of fourth order method to numerically solve y'=f(x,y) for initial x0,y0 going to x1 with n increments, returns y at x1")
function RungeKutta(f,x0,y0,x1,n) = (
	# Note we can't check the 2 arguments, FIXME
	if not IsFunction(f) then
		(error("RungeKutta: f must be a function of two arguments");bailout)
	else if not IsValue(x0) or not IsValue(y0) or not IsValue(x1) then
		(error("RungeKutta: x0, y0 and x1 must be numbers");bailout)
	else if not IsPositiveInteger(n) then
		(error("RungeKutta: n must be a positive integer");bailout);
	h := float(x1-x0) / n;
	x := float(x0);
	y := float(y0);
	for k = 1 to n do (
		# This is unoptimized
		k1 := f(x,y);
		x1 := x + h/2;
		y1 := y + k1*h/2;
		k2 := f(x1,y1);
		x2 := x1;
		y2 := y + k2*h/2;
		k3 := f(x2,y2);
		x3 := x + h;
		y3 := y + k3*h;
		k4 := f(x3,y3);
		y := y + (1/6)*(k1 + 2*k2 + 2*k3 + k4)*h;
		x := x + h
	);
	y
)
protect ("RungeKutta")