System of three coupled linear differential equations Edit: Added a new part to the bottom after doing some thinking
Good day, 
I am trying to solve the following system of equations, where all the M's depend on time.
$\dot{M_1} = -p_1M_1 + p_2M_2$
$\dot{M_2} = -(p_2+p_3)M_2 + p_1M_1+p_4M_3$
$\dot{M_3} = -p_4M_3 + p_3M_2$
Now to be honest, I am rather clueless as for where to start. I have solved such a system once before, but that was using an adiabatic approximation, e.g. saying that one of the differential equations was approximately zero on the timescale at which the others change. I don't really have such information now. 
However, from the context from which I get the problem I can say that at a certain time, I am in one of the three M's and not in the others. I suppose I could go from there, as a sort of initial condition, and try to go further? 
Or perhaps it would be possible to make a simulation of the dynamics instead. I have both matlab and mathematica available, but I an inexperienced with both. This does not mean that I am not able to read documentation of course, so I'll have to look into that. I don't really know where to start, but it seems very plausible that this can be simulated.
After working on it for a while, with the help of the other posters
Alright, so I'm not sure if the answer is correct. My answer comes out to be 


But this can't be right, can it? Because here it seems that if I start in the situation where M1(0) = 1, M2(0)=M3(0) = 0, then the dynamics are very simple: M3(t) is just equal to M1(0) for all time, M1(t) simply exponentially decays, and M2(t) exponentially grows. 
This is in no way coupled at all, while they are supposed to depend on each other all throughout their time evolution.
I wrote a matlab script that solves the system for me, and with that it does show the correct dynamics. With equal probabilities, it converges to all three being 0.333 at high times, which is what I want. However, the solver doesn't give me the actual formula's, so I don't see what I did wrong.
My code:
%differential numerically
p1 = 0.01;
p2 = 0.05;
p3 = 0.05;
p4 = 0.05;
syms M1(t) M2(t) M3(t)
Y = dsolve(diff(M1) == -p1*M1+p2*M2, diff(M2) == -(p2+p3)*M2+p1*M1+p4*M3, diff(M3) == -p4*M3+p3*M2, M1(0) == 1, M2(0) == 0, M3(0) ==0);
t=linspace(0,500,20000);
M1=eval(vectorize(Y.M1));
M2=eval(vectorize(Y.M2));
M3=eval(vectorize(Y.M3));
figure
plot(t,M1,'b-',t,M2,'g-',t,M3,'r-');

 A: You might consider the matrix form of the problem, since it is linear. Indeed, if $\vec{M} = M_i$, then:
$$\frac{d\vec{M}}{dt} = A \vec{M} =  \left( \begin{array}{ccc} -p1 & p_2 & 0 \\ p_1 & -(p_2+p_3) & p_4 \\ 0 & p_3 & -p_4 \end{array}  \right) \, \vec{M}.$$
Hence, the matrix $\Phi(t) = e^{t A}$ (exponential matrix) is said to be a fundamental matrix of the sistem, and the solution becomes:
$$M(t) = \Phi(t) \, \vec{C},$$
where $\vec{C} = C_i$ is the vector of constants of integration.
Cheers.
A: You say 

This is in no way coupled at all, while they are supposed to depend on each other all throughout their time evolution.

But this is alright - you hand-picked a particular initial condition that made one of the quantities constant. In general though, the system will still be coupled.
This is related to the fact that, for many linear operators, there is a choice of basis so that the matrix of the linear operator is diagonal.
For instance, given the system
$$
\dot x = 2 x + y, \\
\dot y = x + 2y,
$$
in matrix form this is
$$
\begin{pmatrix}\dot x \\ \dot y\end{pmatrix} = \begin{pmatrix}2 & 1 \\ 1 & 2\end{pmatrix}\begin{pmatrix}x \\ y\end{pmatrix}
$$
which has general solution
$$
\begin{pmatrix}x(t) \\ y(t)\end{pmatrix} = c_1\begin{pmatrix}1 \\ -1\end{pmatrix}e^t + c_2 \begin{pmatrix}1 \\ 1\end{pmatrix}e^{3t}
$$
and hopefully there's no argument that $x$ and $y$ still look coupled.
However, if we first change variables to 
$$
w = x + y, \\
z = x - y,
$$
then the system can be written
$$
\dot w = \dot x + \dot y = 3x + 3y = 3w \\
\dot z = \dot x - \dot y = x - y = z
$$
and we have effectively decoupled the system; i.e. the corresponding matrix equation would be diagonal, and it's pretty clear that $w$ and $z$ are independent.
EDIT: to compute the matrix exponential of your system (let's assume all the $p_i$'s are equal), find the eigenvalues and eigenvectors of 
$$
A = \begin{bmatrix}-p & p & 0 \\ p & -2p & p \\ 0 & p & -p \end{bmatrix} = p\begin{bmatrix}-1 & 1 & 0 \\ 1 & -2 & 1 \\ 0 & 1 & -1 \end{bmatrix}
$$
which are $0, -p, -3p$. The corresponding eigenvectors are
$$
\begin{pmatrix}1 \\ 1 \\ 1\end{pmatrix}, \begin{pmatrix}1 \\ 0 \\ -1\end{pmatrix}, \begin{pmatrix}1 \\ -2 \\ 1\end{pmatrix},
$$
so solutions are 
$$
\begin{pmatrix}1 \\ 1 \\ 1\end{pmatrix}, \begin{pmatrix}1 \\ 0 \\ -1\end{pmatrix}e^{-pt}, \begin{pmatrix}1 \\ -2 \\ 1\end{pmatrix}e^{-3pt},
$$
We could get a fundamental matrix by taking these solutions as columns, but we want a fundamental matrix so that at $t=0$ we get the identity. One way to do this is to solve initial value problems for $(1,0,0)^T$, $(0,1,0)^T$ and $(0,0,1)^T$, then use those solutions as the columns.
