\ Numerical solution of first-order ordinary differential equation
\ by 4th order Runge-Kutta algorithm following
\ Abramowitz & Stegun p. 896, 25.5.10


\ ---------------------------------------------------
\     (c) Copyright 1999  Julian V. Noble.          \
\       Permission is granted by the author to      \
\       use this software for any application pro-  \
\       vided this copyright notice is preserved.   \
\ ---------------------------------------------------

\ This is an ANS Forth program requiring the
\   FLOAT, FLOAT EXT, FILE and TOOLS EXT wordsets.

\ Environmental dependencies:
\       Assumes independent floating point stack


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

\    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  ;
9 fvariables t h h2 x tf  xk1 xk2 xk3 xk4


v: fdummy

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


: 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
        f" h2 = h/2"
        BEGIN   display   x'        \ indefinite loop
        done?   UNTIL  ;