\begin{align} \dot E &= f_E(E,L,A) = -\mu_e E - c_{le} E L - c_{ae} E A - a_e E + b A \\ \dot L &= f_L(E, L) = -\mu_L L + a_e E - a_L L \\ \dot P &= f_P(L,P) = -\mu_P P - a_p P + a_L L\\ \dot A &= f_A(P,A) = -\mu_A A + a_p P \end{align}
Has the corresponding variance-covariance dynamics
Where g_i is the second jump moment, which is a function of the state (E, L, P, A) just as f_i is. (In this case it will correspond to the sum of all birth and death terms).
General Form & Algorithm
Consider the dynamics are given by
\dot X_i = f_i(\vec X)
and define variance-covariance matrix M and the Jacobian matrix of f as J. Then the dynamics of the diagonal elements (variance terms) are written as