Software to find the equilibria of a set of differential equations? Context for this question: I'm in mathematical biology, and while I've got a decent grasp of the concepts, and am a decent hand with simulation and numerical approaches, my analytical abilities and general skill at manipulating equations has always been rough. Probably the root of deciding I was "bad at math" back in the day, though that's neither here nor there.
What I'm trying to do is find the equilibriums for a set of differential equations, similar to the following:
$$\begin{align}
\frac{{dS}}{{dt\strut}} =& \mu  - (\beta I + \mu )S\\
\frac{{dE}}{{dt\strut}} =& \beta SI - (\mu  + \sigma )E\\
\frac{{dI}}{{dt\strut}} =& \sigma E - (\mu  + \gamma )I\\
\frac{{dR}}{{dt\strut}} =& \gamma I - \mu R
\end{align}$$
I can do this one, though admittedly it takes me quite some time, but more complex systems involve a fair number of hours going down dead ends that really shouldn't be, basic errors etc. 
Any chance any of the myriad math software packages out there have a way to simplify the computation/manipulation steps of finding the equilibria for a system like this? Any general tips if there aren't? You've got Matlab, Mathematica, and really any open source package you please to work with.
 A: Since at the equilibrium time derivatives are zero, the issue reduces to solving a system of algebraic equations. Mathematica has many commands for equation solving, among them Solve and Reduce.
Here is a snapshot of solving your particular system:

A: In Maple, it's pretty straightforward.

eqs:= {mu - (beta*i + mu)*S = 0,
         beta*S*i - (mu+sigma)*E = 0,
         sigma*E - (mu+gamma)*i = 0,
         gamma*i - mu*R = 0};
  solve(eqs,{S,E,i,R});

{E = 0, R = 0, S = 1, i = 0}, {E = -mu*(-beta*sigma+mu^2+mu*gamma+sigma*mu+sigma*gamma)/sigma/(mu+sigma)/beta, R = -gamma*(-beta*sigma+mu^2+mu*gamma+sigma*mu+sigma*gamma)/beta/(mu^2+mu*gamma+sigma*mu+sigma*gamma), S = (mu^2+mu*gamma+sigma*mu+sigma*gamma)/beta/sigma, i = -mu*(-beta*sigma+mu^2+mu*gamma+sigma*mu+sigma*gamma)/beta/(mu^2+mu*gamma+sigma*mu+sigma*gamma)}
A: To find Equilibrium or Steady State of Differential Equations using MATLAB would be,
    clc, clear all;
    syms S E I R mu beta sigma gamma 'real'
    e1 = eval('mu-(beta*I+mu)*S');
    e2 = eval('beta*S*I-(mu+sigma)*E');
    e3 = eval('sigma*E-(mu+gamma)*I');
    e4 = eval('gamma*I-mu*R');
    sol=solve(e1,S,e2,E,e3,I,e4,R);
    disp([sol.S sol.E sol.I sol.R])

Output;
$$
S=(\gamma\mu + \gamma\sigma + \mu\sigma + \mu^2)/(\beta\sigma),\ \\ E=-(\gamma\mu^2 + \mu^2\sigma +\mu^3 - \beta\mu\sigma + \gamma\mu\sigma)/(\beta\sigma^2 + \beta\mu\sigma),\ \\ I=-(\gamma\mu^2 + \mu^2\sigma + \mu^3 - \beta\mu\sigma + \gamma\mu\sigma)/(\beta\mu^2 + \beta\gamma\mu + \beta\gamma\sigma + \beta\mu\sigma),\ \\ R=-(\gamma^2\mu + \sigma\gamma^2 + \gamma*\mu^2 + \sigma\gamma\mu - \beta\sigma\gamma)/(\beta\mu^2 + \beta\gamma\mu + \beta\gamma\sigma + \beta\mu\sigma)
$$ 
