Software for numerical solution of a non-linear ODE system? 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.
 A: 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.
A: Mimicking user110985 s code in octave, first define a function in it's own file:
function lRet = xdot_se1 (x, t)
  lambda = 1;
  k = 1;
  n = 3;
  xdot    = zeros(3,1);
  xdot(1) = x(2);
  xdot(2) = -lambda/(x(3).^(n) - k^2*lambda + eps);
  xdot(3) = 2*n*x(1)*x(2)*lambda*x(3)^n ./ (eps + 2*n*x(3)^(2*n-1) - 2*n*lambda*k^2*x(3)^(n-1) -n*(n-1)*x(1)^2*lambda*x(3)^(n-1));
  lRet = xdot;
endfunction`

Then a script file doing the following:
close all; clear all; clc;

t = (0:0.05:5)';
a = lsode(@xdot_se1,[1,1,4]',t');

figure(1);

plot(t',a(:,1),'r','linewidth',2); hold on;
plot(t',a(:,2),'g','linewidth',2);
plot(t',a(:,3),'b','linewidth',2);`

We get the similar plot:

