Zane (and everyone else), You might be interested in trying VFGEN (www.warrenweckesser.net/vfgen). I just released version 2.2.0. VFGEN is a program that takes an XML description of a vector field (i.e. differential equations) and generates code for a wide variety of numerical software packages, including SciPy. A simple example is the damped pendulum: θ' = v v' = -b*v/(m*L^2) - (g/L)*sin(θ) An XML vector field file for this system is ---------------------------------------------------------------- <?xml version="1.0" ?> <VectorField Name="pendulum" Description="Pendulum Vector Field"> <Parameter Name="g" Description="gravitational constant" DefaultValue="9.81" /> <Parameter Name="b" Description="friction constant" DefaultValue="0.0" /> <Parameter Name="L" Description="pendulum length" DefaultValue="1.0" /> <Parameter Name="m" Description="mass" DefaultValue="1.0" /> <StateVariable Name="theta" Description="Angle, measured from straight down" Formula="v" PeriodFrom="0" PeriodTo="2*Pi" DefaultInitialCondition="Pi-0.01" /> <StateVariable Name="v" Description="angular velocity" Formula="-b*v/(m*L^2)-(g/L)*sin(theta)" DefaultInitialCondition="0.0" /> <Function Name="energy" Description="total energy (kinetic plus potential)" Formula="m*L^2*v^2/2 - m*g*L*cos(theta)" /> </VectorField> ---------------------------------------------------------------- The VFGEN scipy command reads this file and generates this Python code: ---------------------------------------------------------------- # # pendulum.py # # Python file for the vector field named: pendulum # This file implements the vector field as # the function pendulum_vf. This function # can be used with the SciPy ODEINT function. # # This file was generated by the program VFGEN (Version:2.2.0) # Generated on 17-Mar-2008 at 00:48 # from math import * import numpy # # The vector field. # def pendulum_vf(y_,t_,p_): """ The vector field function for the vector field "pendulum" """ Pi = pi theta = y_[0] v = y_[1] g = p_[0] b = p_[1] L = p_[2] m = p_[3] f_ = numpy.zeros((2,)) f_[0] = v f_[1] = -sin(theta)/L*g-1.0/m/(L*L)*b*v return f_ # # The Jacobian. # def pendulum_jac(y_, t_, p_): """ The Jacobian of the vector field "pendulum" """ Pi = pi theta = y_[0] v = y_[1] g = p_[0] b = p_[1] L = p_[2] m = p_[3] # Create the Jacobian matrix: jac_ = numpy.zeros((2,2)) jac_[0,1] = 1.0 jac_[1,0] = -1.0/L*cos(theta)*g jac_[1,1] = -1.0/m/(L*L)*b return jac_ # # User function: energy # def pendulum_energy(y_, t_, p_): """ The user-defined function "energy" for the vector field "pendulum" """ Pi = pi theta = y_[0] v = y_[1] g = p_[0] b = p_[1] L = p_[2] m = p_[3] return m*(L*L)*(v*v)/2.0-m*L*cos(theta)*g ---------------------------------------------------------------- VFGEN can also generate a simple command-line solver written in Python that uses the above function to solve the initial value problem. Take a look here: http://www.warrenweckesser.net/vfgen/menu_scipy.html If you decide you would like to use a different package (e.g. MATLAB, Scilab, C with CVODE or with the GNU Scientific Library, Fortran with LSODA or RADAU5, etc), another one-line VFGEN command will generate the code that you need. Best regards, Warren Weckesser