I have a 9 systems of nonlinear ODEs to solve. I want to determine the endemic equilibrium points.

How do I go about it?

I tried manual calculation but it becomes cumbersome. Can software be used for this purpose?

  • $\begingroup$ Do you need an analytical solution or will it suffice to have a numerical solution for a particular choice of parameters in your ODE? $\endgroup$
    – Paul
    Aug 7, 2017 at 18:29
  • $\begingroup$ Yes, I will prefer an analytical solution.. N which can you help with? $\endgroup$ Aug 7, 2017 at 18:34
  • $\begingroup$ With analytical solutions, I can then see the behaviour of the system as parameters are being altered.. $\endgroup$ Aug 7, 2017 at 18:36
  • $\begingroup$ Use Mathematica or something like SymEngine in Julia. There's a good change there is no analytical solution though, so numerical solutions may be the only option. $\endgroup$ Aug 7, 2017 at 20:24
  • $\begingroup$ Endemic equilibrium(i.stack.imgur.com/DQq8v.jpg Here s d equation)](i.stack.imgur.com/DQq8v.jpg Here s d equation) $\endgroup$ Aug 8, 2017 at 11:12

1 Answer 1


If you want to find the equilibrium points you need to write your system as a first-order one

$$\frac{d\mathbf{x}}{dt} = \mathbf{F}(\mathbf{x})$$

and solve the non-linear system of equations

$$\mathbf{F}(\mathbf{x})=\mathbf{0}\ ,$$

that might not be solvable analytically. After that, you need to find the eigenvalues of the Jacobian matrix around each equilibrium point $\mathbf{x}^*$, $\mathbf{J}(\mathbf{x}^*)$. You can use a Computer Algebra System (CAS) for the job. SymPy, Maxima, Mathematica and Maple are some examples.

An example of this, using SymPy, for the Lotka-Volterra model is the following

from sympy import symbols, Matrix, solve, diff
x, y = symbols("x y")
alpha, beta, gamma, delta = symbols("alpha beta gamma delta")
dXdt = Matrix([alpha*x - beta*x*y, delta*x*y -gamma*y])
eqpts = solve(dXdt, [x, y])
J = Matrix(2, 2, lambda i, j: diff(dXdt[i], [x, y][j]))
vals1 = J.subs({x: eqpts[0][0], y: eqpts[0][1]}).eigenvals()
vals2 = J.subs({x: eqpts[1][0], y: eqpts[1][1]}).eigenvals()

with output

{-gamma: 1, alpha: 1}
{sqrt(-alpha*gamma): 1, -sqrt(-alpha*gamma): 1}

Here you should take into account two things:

  1. The system of non-linear equations might not have an analytical solution, and you might need to solve it numerically (with Newton method, for example).
  2. Even if you have the analytical solution(s) for your system, the eigensystem might not be solvable analytically.

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.