\ Numerical solution of first-order ordinary differential equation
\ by 2nd order Runge-Kutta algorithm following
\ Abramowitz & Stegun p. 896, 25.5.6

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

\    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 )runge





MARKER  -runge

INCLUDE ftran111.f


: fvariables     0 DO  FVARIABLE  LOOP  ;
5 fvariables t h x tf  xk


v: fdummy

: x'    \ integration step
        f" xk = h*fdummy(x,t)"          \ compute k
        f" t = t+h"                     \ increment t
        f" xk + h*fdummy(x+xk,t)"  F2/  ( f: -- dx)
        x F@   F+   x F!                \ x = x + dx
;


: 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>  ;


: )runge   ( xt --)  ( f: x0 t0 tf h -- )
        defines fdummy              \ vector fn_name
        h F!  tf F!  t F!  x F!     \ initialize variables
        BEGIN   display   x'        \ indefinite loop
        done?   UNTIL  ;