package require math::bigfloat

 namespace import ::math::bigfloat::*

 # compute y'=f(t,y) with the initial value y(0)
 # here t0..tn and n are the steps and their number
 # f is the function (with two arguments : t and y)
 # y0 is y(0)
 proc RungeKutta4 {t0 tn f y0 n} {
    # f(x_0)=current value
    set t $t0
    # y(initial_value)
    set y $y0
    set step [div [sub $tn $t0] [fromstr $n]]
    set half_step [div $step [fromstr 2]]
    # from x_1 to x_N
    for {set i 0} {$i < $n} {incr i} {
        set k1 [$f $t $y]
        set k2 [$f [add $t $half_step] [add $y [mul $half_step $k1]]]
        set k3 [$f [add $t $half_step] [add $y [mul $half_step $k2]]]
        set k4 [$f [add $t $step] [add $y [mul $step $k3]]]
        set slope [div [add [div [add $k1 $k4] [fromstr 2]] [add $k2 $k3]] [fromstr 3]]
        set y [add $y [mul $slope $step]]
        set t [add $t $step]
    }
    # the return value is a bigfloat
    return $y
 }

 proc ident {t y} {
    set y
 }

 # y'(t)=y(t) and y(0)=1

 set x0 [fromstr 0.0 10]
 set y0 [fromstr 1.0 10]
 set xEnd [fromstr 1.0 10]
 set steps 100
 puts "With RK4 exp(1)=[tostr [RungeKutta4 $x0 $xEnd ident $y0 $steps]]"
 puts "Real value equals : [tostr [exp $xEnd]]"