\ Numerical solution of first-order ordinary differential equation
\ by 4th order Milne predictor-corrector algorithm following
\ Abramowitz & Stegun p. 896, 25.5.13

\ Solves      dx/dt = f(x,t)

\ P:    x(n+1) = x(n-3) + (4h/3) * [2f(n) - f(n-1) +2f(n-2)]
\ C:    x(n+1) = x(n-1) + (h/3)  * [f(n-1) + 4f(n) + f(n+1)]

\    This is an ANS Forth program requiring FLOATING and FLOATING EXT
\    wordsets. This program assumes a separate floating point stack.

\    For simplicity and clarity, the program is intended to be compiled
\    using a FORmula TRANslator


\ Example
\    : fnb   ( f: x t -- t^2*exp[-x])
\        f^2  FSWAP  FNEGATE  FEXP  F*  ;
\    use( fnb 0e0 0e0 5e0 0.1e0 )milne



MARKER  -milne

INCLUDE ftran111.f


: fvariables     0 DO  FVARIABLE  LOOP  ;

5 fvariables t h h3 h43 tf
8 fvariables x x0 xm1 xm2 xm3 y y0 ym1
3 fvariables f0 fm1 fm2
4 fvariables xk1 xk2 xk3 xk4


v: fdummy

: update    \ shuffle variables
    f" xm3 = xm2"
    f" xm2 = xm1"                   f" fm2 = fm1"
    f" xm1 = x0"    f" ym1 = y0"    f" fm1 = f0"
    f" x0 = x"      f" y0 = y"
;


: runge_step
        f" f0 = fdummy(x,t)"
        f" xk1 = h*f0/2 "                  \ compute k1
        f" t = t  +h/2 "                   \ t -> t + h/2
        f" xk2 = h*fdummy(x+xk1,t)"
        f" xk3 = h*fdummy(x+xk2/2,t)"
        f" t = t+h/2 "                     \ t -> t + h/2
        f" xk4 = h*fdummy(x+xk3,t)/2"
        f" x = x + (xk1 + xk2 + xk3 + xk4)/3 "
;

: x'    \ integration step
    f" f0 = fdummy(x,t)"
    f" x = xm3 + h43 * (2*(f0 + fm2) - fm1)"
    f" t = t + h"
    f" y = ym1 + h3 * (fm1 + 4* f0 + fdummy(x,t))"
;


: display       CR   t F@  FS.  5 SPACES   x F@  FS.
                                5 SPACES  f" ln(1+t^3/3)"  FS.
;


: done?  t F@  tf F@  F>  ;


: initialize    ( xt --)  ( f: x0 t0 tf h -- )
        defines fdummy              \ vector fn_name
        h F!  tf F!  t F!  x F!     \ initialize variables
        f" h3 = h/3"
        f" h43 = 4*h3"
        0.e0  FDUP  xm2 F!          \ put 0 in these
              FDUP  xm1 F!
              FDUP  x0  F!
              FDUP  y0  F!
                    fm1 F!
        4 0 DO   runge_step     f" y = x"   update   LOOP
;

: )milne   ( xt --)  ( f: x0 t0 tf h -- )
        initialize
        BEGIN   display   x'        \ indefinite loop
        done?   UNTIL  ;