[Arjen Markus] (5 march 2009) Yesterday I spoke with [Kevin Kenny] and [Cameron Laird] about the possibilities
of developing an application a la XPPAUT [http://www.math.pitt.edu/~bard/xpp/xpp.html]. The idea is to
use such packages as [Tclode] and [Plotchart] to solve systems of ordinary differential equations and plot 
the results.

These equations should be definable in a very similar way as XPPAUT allows. The user-interface should be
much more modern though.

Anyway, here is a (second) shot at something like that. I have called it "Bernoulli" in honour of those 
great 17th-century mathematicians, Johan, Jacob and Daniel. (It is also a bit of nostalgia for me, as 
some 15-20 years ago I developed a similar program, also called bernoulli).

The version below defines a number of commands and can actually integrate systems of ODEs using the
Tclode extension. No neat user-interface yet, just the basics. See the example at the end for how to
use the commands. (The plotnow command simply plots one variable as function of the other after a 
computation has been done. Each time a new plot is created)

----
======
# bernoulli.tcl --
#     Integrate systems of ordinary differential equations
#
package require Tk
package require Plotchart
if {0} {
set dir tclode1.0
source $dir/pkgIndex.tcl
}
package require Tclode

# namespaces --
#     Create a number of private namespaces
#
namespace eval ::time   {
    variable steps 100
    variable tmin  0.0
    variable tmax  1.0
    variable statevars {t}
}
namespace eval ::series {}
namespace eval ::init   {}
namespace eval ::deriv  {}

# plotf --
#     Convenience procedure to plot a function of one variable
#
# Arguments:
#     range       Range over which to plot the function (x-axis)
#     expression  Expression by which the function is defined
#                 (independent variable: t)
#
# Result:
#     None
#
proc plotf {range expression} {

    # Very simple for now: no parameters

    foreach {tmin tmax} $range {break}
    set nosteps 100
    set tstep [expr {($tmax-$tmin)/double($nosteps)}]

    if { ![winfo exists .c] } {
        pack [canvas .c]
    } else {
        .c delete all
    }

    set tvalues {}
    set xvalues {}
    set xmin    {}
    set xmax    {}

    for {set i 0} {$i < $nosteps+1} {incr i} {
        set t [expr {$tmin + $i*$tstep}]
        set x [expr $expression]
        lappend tvalues $t
        lappend xvalues $x

        if { $xmin == {} || $x < $xmin } {
            set xmin $x
        }
        if { $xmax == {} || $x > $xmax } {
            set xmax $x
        }
    }

    set tscale [::Plotchart::determineScale $tmin $tmax]
    set xscale [::Plotchart::determineScale $xmin $xmax]

    set p [::Plotchart::createXYPlot .c $tscale $xscale]

    foreach t $tvalues x $xvalues {
        $p plot function $t $x
    }
}

# plotnow --
#     Convenience procedure to plot one variable as function of the other
#
# Arguments:
#     xvar        Variable on the x-axis
#     yvar        Variable on the y-axis
#
# Result:
#     None
#
proc plotnow {xvar yvar} {

    if { ![winfo exists .c] } {
        pack [canvas .c]
    } else {
        .c delete all
    }

    set xvalues [set ::series::$xvar]
    set yvalues [set ::series::$yvar]
    set xmin    {}
    set xmax    {}
    set ymin    {}
    set ymax    {}

    foreach x $xvalues y $yvalues {

        if { $xmin == {} || $x < $xmin } {
            set xmin $x
        }
        if { $xmax == {} || $x > $xmax } {
            set xmax $x
        }

        if { $ymin == {} || $y < $ymin } {
            set ymin $y
        }
        if { $ymax == {} || $y > $ymax } {
            set ymax $y
        }
    }

    set xscale [::Plotchart::determineScale $xmin $xmax]
    set yscale [::Plotchart::determineScale $ymin $ymax]

    set p [::Plotchart::createXYPlot .c $xscale $yscale]

    foreach x $xvalues y $yvalues {
        $p plot function $x $y
    }
}

# timeframe --
#     Define the time frame for the next computation
#
# Arguments:
#     tmin        Minimum time
#     tmax        Maximum time
#
# Result:
#     None
#
proc timeframe {tmin tmax} {

    set time::tmin $tmin
    set time::tmax $tmax
}

# init --
#     Define the initial value of a variable
#
# Arguments:
#     varname     Name of the (state) variable
#     value       Initial value
#
# Result:
#     None
#
proc init {varname value} {

    set init::$varname $value
}

# deriv --
#     Register the expression that determines the derivative
#
# Arguments:
#     varname     Name of the (state) variable
#     expression  Expression to be used
#
# Result:
#     None
#
proc deriv {varname expression} {

    if { [lsearch $::time::statevars $varname] < 0 } {
        lappend ::time::statevars $varname
    }
    set deriv::$varname $expression
}

# go --
#     Do the computation
#
# Arguments:
#     None
#
# Result:
#     None
#
proc go {} {

    #
    # Time administration and initial values
    #
    set t           $::time::tmin
    set ::time::dt  [expr {($::time::tmax-$t)/double($::time::steps)}]

    set args {}
    foreach v [lrange $::time::statevars 1 end] {
        lappend args $v [set ::deriv::$v]
        set $v [set ::init::$v]
    }

    #
    # Construct the ODE solver
    #
    set ::time::solver [::tclode::odesolv -atol 1.0e-6 -rtol 1.0e-6 $args]

    #
    # Run the computation
    #
    $::time::solver run $t

    for {set ::time::i 0} {$::time::i <= $::time::steps} {incr ::time::i} {
        set ::time::tnext [expr {$::time::tmin + $::time::dt * $::time::i}]
        $::time::solver continue $::time::tnext
        foreach v $::time::statevars {
            lappend ::series::$v [set $v]
        }
    }
}

# main --
#     Test the stuff
#

#plotf {0 10} {exp(-$t)*cos(2.0*$t)}
#after 3000 {
#    plotf {0.0001 100} {sin($t)/$t} ;# Avoid t=0 of course
#}

init s 0.0
init c 1.0
init x 0.0
deriv s {$c}
deriv c {-$s}
deriv x {1}
timeframe 0.0 10.0
go

catch {
console show
}
foreach t $::series::t s $::series::s c $::series::c x $::series::x {
    puts [format "%10.4f %10.4f %10.4f %10.4f" $t $s $c $x]
}

plotnow t c 
======
----
!!!!!!
%| [Category Mathematics] |%
!!!!!!