Let's see an example of a first order ODE system : $$\begin{align*} y_1'&=a_{11}\cdot y_1+a_{12}\cdot y_2 \qquad(1)\\ y_2'&=a_{21}\cdot y_1+a_{22}\cdot y_2 \qquad(2) \end{align*}$$ where $y_1$ and $y_2$ are functions we are to solve, and $a_{11},\; a_{12},\; a_{21},\; a_{22}$ are constants. I know 2 methods to solve this, one is to eliminate $y_1$ or $y_2$ by replacing, and then we will get 2 decoupled second order ODEs

$$ y_1''-(a_{11}+a_{22})y_1'+(a_{11}a_{22}-a_{12}a_{21})y_1=0 $$ the other ODE is similar.

The other method: Try to calculate the eigenvalues and eigenvectors of the coefficient matrix $\{\{a_{11},a_{12}\},\{a_{21},a_{22}\}\} $, and then get the general solutions for system of $y_1$ and $y_2$ simultaneously. Though the second method is more popular in the textbooks, I like the first one more.

So my question is: Which function in Mathematica can decouple the ODE system by the first method automatically?

  • $\begingroup$ In your 2nd equation isn't the first y2 supposed to be a y1 or is it correct as is? $\endgroup$ Nov 4, 2012 at 13:30
  • $\begingroup$ Changed to y1. Please roll back if I guessed incorrectly. $\endgroup$ Nov 4, 2012 at 13:36
  • $\begingroup$ It seems to me that going to the eigen-space precisely decouples the equations, so the two methods are equivalent(?) $\endgroup$
    – chris
    Nov 4, 2012 at 14:35
  • $\begingroup$ Do I understand correctly that this question is about manipulating the ODE's for further manual analysis? You don't want to solve them with the help of Mathematica, right? $\endgroup$
    – sebhofer
    Nov 4, 2012 at 18:06

2 Answers 2


I believe going to the eigen-space precisely decouples the equations, so the two methods are strictly equivalent.

Lets define your matrix

 mat = {{a11,a12},{a21,a22}};

Now find the eigenvector matrix P (the so called transformation matrix)

 {\[Lambda]s, P} = Eigensystem[mat]; P = Transpose[P];

Let us check that going to the eigen-space makes the Matrix diagonal

 pmat= Inverse[P].mat.P // FullSimplify

 1/2 (a11+a22-Sqrt[(a11-a22)^2+4 a12 a21])  0
 0  1/2 (a11+a22+Sqrt[(a11-a22)^2+4 a12 a21])

Check that the diagonal elements are the eigenvalues

pmat[[1, 1]]/\[Lambda]s[[1]] // Simplify
pmat[[2, 2]]/\[Lambda]s[[2]] // Simplify

So if we define as new variables {z1[t],z2[t]} so that

Thread[{z1[t],z2[t]}=pmat.{y1[t],y2[t]} ]

The system in z is simply

z1'[t] ==  \[Lambda]s[[1]] z1[t]
z2'[t] ==  \[Lambda]s[[1]] z2[t]

which is decoupled.


Not really sure why you would want to prescribe the method Mathematica uses to solve the coupled ODEs. If you don't Mathematica solves them just fine:

eq1 = y1'[t] == a11*y1[t] + a12*y2[t];
eq2 = y2'[t] == a21*y1[t] + a22*y2[t];

FullSimplify[DSolve[{eq1, eq2}, {y1[t], y2[t]}, t]]

Mathematica graphics

But if you want to go the eliminate & replace route you could do the following:

eq1D = Assuming[{a11 != 0, a12 != 0, a21 != 0, a22 != 0}, 
         FullSimplify[eq2 /. Solve[eq1, y2[t]] /. D[Solve[eq1, y2[t]], t]]][[1, 1]]

Mathematica graphics

eq2D = Assuming[{a11 != 0, a12 != 0, a21 != 0, a22 != 0}, 
         FullSimplify[eq1 /. Solve[eq2, y1[t]] /. D[Solve[eq2, y1[t]], t]]][[1, 1]]

Mathematica graphics

DSolve[eq1D, y1[t], t]

Mathematica graphics

DSolve[eq2D, y2[t], t]

Mathematica graphics


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.