It's generally nice and friendly to show the complete code:
In [18]: cP1, cE1, cE2, cE3 = symbols('cP1, cE1:4', cls=Function)
In [19]: t, k1, k6, k8 = symbols('t, k1, k6, k8')
In [20]: eqs = [Eq(cP1(t).diff(t), k1*cE1(t)**3), Eq(cE1(t).diff(t), -k1 * cE1(t)**3 + k6 * cE3(t)**2),
...: Eq(cE2(t).diff(t), -k8 * cE2(t)), Eq(cE3(t).diff(t), k8 * cE2(t) - k6 * cE3(t)**2)]
...:
In [21]: for eq in eqs:
...: pprint(eq)
...:
d 3
──(cP₁(t)) = k₁⋅cE₁ (t)
dt
d 3 2
──(cE₁(t)) = - k₁⋅cE₁ (t) + k₆⋅cE₃ (t)
dt
d
──(cE₂(t)) = -k₈⋅cE₂(t)
dt
d 2
──(cE₃(t)) = - k₆⋅cE₃ (t) + k₈⋅cE₂(t)
dt
In [22]: dsolve(eqs)
---------------------------------------------------------------------------
NotImplementedError Traceback (most recent call last)
<ipython-input-22-69ab769a7261> in <module>
----> 1 dsolve(eqs)
~/current/sympy/sympy/sympy/solvers/ode/ode.py in dsolve(eq, func, hint, simplify, ics, xi, eta, x0, n, **kwargs)
609 "number of functions being equal to number of equations")
610 if match['type_of_equation'] is None:
--> 611 raise NotImplementedError
612 else:
613 if match['is_linear'] == True:
NotImplementedError:
This means that dsolve can not yet handle this particular type of system. Note that in general nonlinear systems of ODEs are very unlikely to have analytic solutions (dsolve is for finding analytic solutions, if you want numerical solutions use something like scipy's odeint).
As nonlinear systems go this one is relatively friendly so it might be possible to solve it. Let's see...
Firstly there is a conserved quantity (the sum of all 4 variables) that we can use to eliminate one equation. Actually that doesn't help as much because the first equation is already isolated from the others: if we knew cE1 we could just integrate but the conserved quantity gives it more easily if the other variables are known.
The structure of the system is like:
cE2 ---> cE3 ---> cE1 ---> cP1
implying that it can be solved as a sequence of ODEs where we solve the 3rd equation for cE2 and then the 4th equation for cE3 and then use that for cE1 and so on. So we can reduce this to a problem involving a sequence of mostly nonlinear single ODEs. That also is very unlikely to have an analytic solution but let's give it a try:
In [24]: dsolve(eqs[2], cE2(t))
Out[24]:
-k₈⋅t
cE₂(t) = C₁⋅ℯ
In [25]: cE2sol = dsolve(eqs[2], cE2(t)).rhs
In [26]: eqs[3].subs(cE2(t), cE2sol)
Out[26]:
d -k₈⋅t 2
──(cE₃(t)) = C₁⋅k₈⋅ℯ - k₆⋅cE₃ (t)
dt
At this point in principle we could solve for cE3 here but sympy doesn't have any way of solving this particular nonlinear ODE so dsolve gives a series solution (I don't think that's what you want) and the only other solver that might handle this is lie_group but that actually fails.
Since we can't get an expression for the solution for cE3 we also can't go on to solve for cE1 and cP1. The ODE that fails there is a Riccati equation but this particular type of Ricatti equation is not yet implemented in dsolve. It looks like Wolfram Alpha gives an answer in terms of Bessel functions:
https://www.wolframalpha.com/input/?i=solve+dx%2Fdt+%3D+e%5E-t+-+x%5E2
Judging from that I guess it's unlikely that we would be able to solve the next equation. At that point Wolfram Alpha also gives up and says "classification: System of nonlinear differential equations":
https://www.wolframalpha.com/input/?i=solve+dx%2Fdt+%3D+e%5E-t+-+x%5E2%2C+dy%2Fdt+%3D+-y%5E3+%2B+x%5E2
I suspect that there is no analytic solution for your system but you could try numerical solutions or otherwise a more qualitative analysis.
