Closed-form solution of a nonlinear system of differential equations I've come up with a nonlinear system of differential equations while analyzing a system approaching chemical equilibrium. This is the equation.
$$
\left\{\begin{aligned}
\frac{dx}{dt}&=px-qx(y+z)-rx(x+y+z)\\
\frac{dy}{dt}&=py-qy(z+x)-ry(x+y+z)\\
\frac{dz}{dt}&=pz-qz(x+y)+r((x+y)^2-z^2)\\
\end{aligned}\right.$$
$x(0)=x_0$, $y(0)=y_0$, $z(0)=0$, $p,q,r$ arbitrary constants
I've tried solving the equations using Wolfram Mathematica and Wolfram|Alpha, yet neither could find a closed-form solution to them. Is it that such solution doesn't exist, or is there a way of finding one, given enough time and computation power?
Btw, I've tried setting $p = 50$, $q = 5$, $r = 5$ and have obtained numerical graphs  of the equations. The curve at the bottom apparently converging to 0 is a plot of $y(t)$, with the one immediately next to it (also converging to 0) being a plot of $x(t)$. The curve that appears to converge to somewhere near 4 is a plot of $z(t)$. However, zooming in a bit more into the curve corresponding to $z(t)$, we find that the function actually is still undergoing oscillation, even for a value of t as high as t = 2.3. 
Thanks!
 A: "Disclosure": When I wrote this, I hadn't seen the differential system.
Observing the curve, and interpretating it as vanishing sinusoidal perturbations superimposed on a decreasing exponential, I obtained the following equation that fits rather well the general form of your curve:
$$\tag{1}y=(3.6+0.4 \cos(2\pi x))e^{-x}$$
(see curve below, starting at $x=0$, not $x=23$).
A more general form will be:
$$\tag{2}y=(A+B \cos(C x))e^{-Dx}$$
with constants that have to be tuned, for example with a "fit" function (like in Matlab). But, at least you have a reasonable track... 
As colleagues that have solved the system dont find solutions that fit the given curve, let us make some retro-engineering in order to try to have hints leading to a rectification of the differential system.
If the solution's shape is like in (2), it means that the eigenvalues of the matrix are $\{-D,-D-iC,-D+iC\}$, with 2 complex eigenvalues, which is not the case here (the eigenvalues of a triangular matrix are its diagonal elements, thus real). Thus, the matrix of the system is, IMHO, not the good one.
A last point: I am puzzled by the scale on the ordinate's axis. Why these repeated $4$'s ?

A: The stationary points with $x^*=0$, $y^*=0$ have from $\frac{dz}{dt}=0$ the third condition  $0=z^*(p−rz^*)$, so that either $z^*=0$ or $z^*=\frac pr=\frac{50}{5}=10$.
The linearization of the system near the fixed point $(x^*,y^*,z^*)=(0,0,p/r)=(0,0,10)$ is
$$
\frac{d}{dt}\pmatrix{x\\y\\z}
=
\pmatrix{
-\frac{pq}r&0&0\\
0&-\frac{pq}r&0\\
-\frac{pq}r&-\frac{pq}r&-p}
\pmatrix{x\\y\\z}
$$
which results in a simple exponential decay with rate $-50$ in all components towards this fixed point. An indeed, the numerical solution confirms that:
import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as pl

p,q,r = 50.0,5.0,5.0;

def func(v,t):
    x,y,z = v
    return [ 
        p*x - q*x*(y+z) - r*x*(x+y+z), 
        p*y - q*y*(z+x) - r*y*(x+y+z),
        p*z - q*z*(x+y) + r*((x+y)**2 - z**2)
        ]    

y0=[ 10.0, 15.0, 0]
x = np.linspace(0,0.25,151)
y=odeint(func, y0, x)

pl.plot(x, y[:,0],"r-."); pl.plot(x, y[:,1],"b-."); pl.plot(x, y[:,2],"g--")
pl.show()


Hypothesis is that you made a typo in implementing the deriatives function.
