# tcl runge-kutta solution of ODES
 # single variable - use rk4step (included for simplicity)
 proc decay { time xn } { ;# simple decay proportional to amount (gives radioactive half life, eg)
  return [expr -2.0*$xn]
 }

 proc rk4step {x time dt dxdtproc } { ;# proc is such that dx/dt = dxdtproc(t)
  set k1 [$dxdtproc  $time  $x] ;# k1 = f (tn, xn)
  set k2 [$dxdtproc  [expr $time+($dt*0.5)] [expr $x + $dt*0.5*$k1]] ;# k2 = f (tn + h?2, xn + h?2 k1)
  set k3  [$dxdtproc [expr $time +$dt*0.5] [expr $x + $dt*0.5*$k2] ] ;# k3= f (tn + h?2, xn + h?2 k2)
  set k4  [$dxdtproc [expr $time +$dt] [expr $x + $dt*$k3] ] ;# k4= f (tn + h, xn + h k3) 
  return  [expr $x + $dt/6.0 *($k1 + 2*$k2 + 2*$k3 + $k4)] ;# h/6(simpsons rule)
 }

 proc test {} {
   set y 0
   set mass 1
   set dt 0.1 ;# try shortening dt to see numerical inaccuracy
 # use single equation R_K scheme - log decay of radioactivity eg (dN/dt=-n-->> N(t)=exp(-t)
   for {set time 0} {[expr $time<=2] } { set time [expr $time+$dt]} {
    puts "At time $time $mass"
    set mass [rk4step $mass $time $dt decay] ;# dmass/dt = -mass - mass=m0(1-exp(time/k))
   }
 }

 test

 # tcl runge-kutta solution of ODES
 # N variables - use rk4step_array. Can be used for higher order ODEs
 # A simple plotting algorithm on a canvas is added.
 set k -0.25
 set omega [expr 2*3.14159]

 # plotted example - damped oscillator;  dv/dt = -k.v -omega.y & dx/dt=v
 proc dxdt { time dt xn} { ;# xn  is array of current variables; dx/dt is first element = speed
        upvar 1 $xn x  ;# array of values
        return $x(0)
 }

 proc oscdvdt { time dt xn} { ;
        upvar 1 $xn x  ;# array of values
        global k omega ;# allows editing, eg to see negatively damped oscilation type "set k 0.2; plot"
        return [expr $k*$x(0)-$x(1)*$omega ] ;# if omega = 2Pi takes 1 sec to oscillate once (Slightly longer if highly damped!)
 }

 proc rk4step_array {vars time dt dxdts } { ;# proc is such that dx/dt = dxdtproc(t,x)
 # for each member of the array 'vars' there is function dxdtproc(i)
 # such that d/dt(vars(i) = dxdtproc(i)
 # dxdtproc(i) has args time, step, array of values
  upvar 1 $vars x  ;# upvar 1 means: use the variable in the calling routine - called $vars with local name x
  upvar 1 $dxdts dxdtproc ;# list of functions for the rhs of equations
  set neqs [array size dxdtproc] ;# number of equations being solved

  for {set i 0} {$i<$neqs} { incr i} {
      set k1($i) [ $dxdtproc($i)  $time $dt x] ;# k1 = f (tn,   xn,   yn)
      set xnu($i) [expr $x($i) + 0.5*$dt*$k1($i) ] ;# next values to evaluate df/dx
  }

 for {set i 0} {$i<$neqs} { incr i} {
         set k2($i) [$dxdtproc($i)   $time $dt xnu]
        set xnu2($i) [expr $x($i) + 0.5*$dt*$k2($i) ] ;# next values
 }

  for {set i 0} {$i<$neqs} { incr i} {
         set k3($i) [$dxdtproc($i) $time $dt xnu2]
        set xnu3($i) [expr $x($i) + $dt*$k3($i) ] ;# next values
  }

  for {set i 0} {$i<$neqs} { incr i} {
           set k4($i) [$dxdtproc($i) $time $dt xnu3]
 }
  for {set i 0} {$i<$neqs} { incr i} {
        set dv [expr  $k1($i) + 2*($k2($i) + $k3($i)) + $k4($i) ]
         set x($i)  [expr $x($i) + $dt*0.1666667 * $dv]
  }
 }

 proc plot {} { ;# uses the simple graph plotter found on https://wiki.tcl-lang.org/8552
  # --------------------------
  #   params
  # --------------------------
  # title
  # canvas width & height
  # delay between plots
  # x = f(t)
  # initial & final times
  array set params \
  {
    title     {Damped Oscillator}
    width     400
    height    400
    delay     0
    x         {$t / 50.}
    t0        0
    t1        400
  }

  # --------------------------
  #   plotting
  # --------------------------
  # computed heights used to scale vertical units
  set h $params(height)
  set h1 [expr {int($h * 0.5)}]  ;# canvas mid-height
  set h2 [expr {$h1 + 1}]
  set h3 [expr {int($h * 0.4)}]  ;# graph mid-height
  # canvas
 if [winfo exists .c] { ;# delete it. Same as clear.
  destroy .c
 }
  canvas .c -width $params(width) -height $h \
            -xscrollincrement 1 -bg beige
  pack .c
  # plotting
  wm title . $params(title)
  set t $params(t0)

 # GWM start up the oscillator solver
  array set fns [list 0 oscdvdt 1 dxdt] ;# dv/dt and dx/dt equations to be solved in 'parallel'.
  array set vars [list 0 0 1 0.5] ;# initial values of time and X

     set tspace [expr $params(t1)/8.0] ;# how often to draw grid
     set nexttic $tspace
      .c create text 20 10 -anchor n -text "Velocity" -fill red
      .c create text 20 25 -anchor n -text "Position" -fill black

  while {$t != $params(t1)} \
  {
    update
    set time [expr $t/30.0]

    rk4step_array vars $time 0.02 fns ;# integrate the second order ODE

    set x [expr $params(x)]
    set vv $vars(1)
    set v [expr {int($vv * $h3) + $h1}]
    if {$t >= $nexttic} {
      set nexttic [expr $nexttic+$tspace]
      .c create text $t 0 -anchor n -text $t -fill gray
      .c create line $t 0 $t $h -fill gray
    }
    .c create line $t $h1 $t $h2 -fill gray
    .c create rectangle $t $v $t $v
    set vv $vars(0)
    set v [expr {int($vv * $h3) + $h1}]
    .c create rectangle $t $v $t $v -outline red
    incr t
  }
 }

 for {set k 0 } {$k>-1.0} {set k [expr $k-0.1]} {
    after 120 ;# wait for a moment....
    plot
 }