[Math] Software for numerical solution of a non-linear ODE system

Question

I have been given a nonlinear system of ODEs which has arisen out of a colleague's engineering research:

$$\begin{array}{rcl}
\dot{x}_0&=&x_1\\
\dot{x}_1&=&-\frac{\lambda}{(x_2)^n-k^2\lambda}x_0\\
\dot{x}_2&=&\frac{2nx_0x_1\lambda(x_2)^n}{2n(x_2)^{2n-1}-2n\lambda k^2(x_2)^{n-1}-n(n-1)(x_0)^2\lambda(x_2)^{n-1}}
\end{array}$$

(I have values for $\lambda$, $k$ and $n$, and also some initial conditions).

Basically all I want is an (open source) computer system to fling them at, and see what happens.

I have so far tried Octave and its "lsode" command (no good; gave errors), Python with "odeint" from sympy.integrate (gave a solution); I need to test out a few others.

I am far from being an expert (or even vaguely competent) at non-linear ODE systems, and I'm hoping for advice as to which system I can use with confidence to generate trustworthy numerical solutions.

Answer 1

I would use Scilab for this job. There are several ODE solvers implemented in Scilab; in the code given below I left the choice of solver to the program, but you can compare the results of different solvers by specifying them in the command.

For numerical values I used $\lambda=k=1$ and $n=3$. Starting point is the column vector $(1,1,4)^T$. Note that Scilab enumerates vector entries starting with $1$, not with $0$ as in your post.

Code:

function xdot=f(t, x)
    xdot=[x(2); -x(1)/(x(3)^3-1); 6*x(1)*x(2)*x(3)^3/(6*x(3)^5-6*x(3)^2-6*x(1)^2*x(3)^2)]
endfunction
x0=[1;1;4]
t0=0
t=0:0.05:5
x = ode(x0,t0,t,f);
plot(t,x(1,:),'r')
plot(t,x(2,:),'g')
plot(t,x(3,:),'b')

Graphical output:

The numerical output is in the matrix x.