DIFFEQ: Solve (Systems of) Differential Equations

convert any higher order differential equations to systems of first order equations ⇾ example, e.g.:
- y1" + a*y1' + k*y1 + F = 0
- convert to
  - y1' = y2
  - y2' = -a*y2 - k*y1 - F
- DIFFEQ(CallBack=subname, T=t, Y=z, DY=dz, T0=t0, Go=go, TOLerance=err, DT1=Step1, itrOK=ok)
Update any f(y,t,...)-dependencies in callback subroutine sub:

required keywords	type	mini sample	Keyword sequence is irrelevant
CallBack	SUB	CB=sub	subroutine sub called back from DIFFEQ
T	NUM	T=x	the independent variable
Y	NUM	Y=y	if only 1 ODE: the dependent variable
	arr	Y=vec	dimension(N) : N ODEs, N dependent variables
DY	arr	DY=dy	same dimension as Y, to be set in callback subroutine sub
optional keywords
Go	num	G=50	max iterations (default is -1: from T0 to T)
T0	num	T0=3	initial value of T (default = 0)
DT1	num	DT1=1e-8	initial stepsize (default = 0.001)
TOLerance	num	TOL=0.01	error limit for each iteration (def = 1E-5)
NonZero	num	NZ=1e-8	y > 0 error control threshold to avoid numeric problems caused by small values of the dependent variable(s)
	0	NZ=0	y error control OFF (allow y < 0 results again)
itrOK	NUM	OK=L	L==1: after a successful iteration, else 0 use for output, or to update any f(y,t,..) work
ERror	LBL	ER=99	on error jump to label 99

solve systems of
parabolic_differential_equations
_or_
reaction-diffusion_equations
- ∂y/∂t = ∂ ²(E⋅y)/∂x ² + ∂(v⋅y)/∂x + R(x, y, y2, ..., t)
Eddy diffusion E, may be a function of t and y, physical dimension e.g. m ³/m/sec, or miles ²/hour
linear velocity v, may be a function of t and y, physical dimension e.g. m ³/m ²/sec, or miles/hour
reaction rate R mole/m ³/sec or cars/mile/hour.
R usually couples multiple reaction-diffusion equations
called without DY ...)
boundary conditions at
- x = xfeed: d(E⋅y)/dx = (v⋅y) - (Vfeed⋅Yfeed)
- x = xout: d(E⋅y)/dx = 0

required keywords	type	mini sample	Keyword sequence is irrelevant
CallBack	SUB	CB=sub	subroutine sub called back from DIFFEQ
T	NUM	T=t	the independent variable. Physical dimension e.g. sec
Y	arr	Y=vec	array of length N: N/M coupled PDEs (N == M if only 1 PDE). Physical dimension e.g. mole/m ³ or cars/mile
X	arr	X=xvec	array of length M: x pillar positions (variable space OK). Physical dimension e.g. m, evaluated at M positions
YFeed	arr	YF=inp	array of length N/M (input at x=0 or x=xmax, see Velocity)
Velocity	arr	V=flow	array of length N (space dependent), N/M (space constant) flow > 0 enters at x=0, flow < 0 at x=xmax
Eddy	arr	E=dif	array of length N (space dependent), N/M (space constant)
Rateatx	arr	R=r	array of length N: differential in/out along X
optional keyword
BCFdiff	num	bcf=0	differences: -1:backward, 0=central, +1=forward