I'm trying to write a Mathematica script to solve this inhomogeneous ODE system of the form $\dot x(t)=M x(t) + g(t)$, where
$M=\left( \begin{array}{ccc} 0 & 0 & 1 \\ 0 & 0 & \frac{1}{l+\frac{1}{2}} \\ -\frac{\text{$\omega $1}^2}{2 \left(\frac{1}{2}-\frac{8 (-2+\pi )}{(2 l+1) \pi ^3}\right)} & -\frac{2}{\left(\frac{1}{2}-\frac{8 (-2+\pi )}{(2 l+1) \pi ^3}\right) \pi } & 0 \\ \end{array} \right) $ and $g(t)=\left(\begin{array}{c}0\\\frac{p(t)}{l+\frac{1}{2}}\\c*p'(t)\end{array}\right)$
The homogeneous part is easy. For the inhomogenous part, my first thought was to use variation of parameters to find a particular solution of the form $M\int M^{-1}g(t)dt$. However this requires $det(M) \neq 0$, which is not satisfied. How wrong would it be to use the pseudoinverse of M in this case? Are there any other methods I can use to find a particular solution of this problem? Eventually I want to scale the script up to a 7-dimensional system, and Mathematica's DSolve function is hanging for larger matrices. Here's the Mathematica code for the matrices if anybody wants to try running it http://pastebin.com/3y7VmZ2v/
It is not the matrix $M$ that is used in the variation of constant formula. Let $Φ(t)$ be a fundamental matrix, i.e., a matrix solution of $$ \dot Φ(t)=M·Φ(t) $$ with $Φ(0)$ non-singular. Then $Φ(t)$ is non-singular everywhere and the particular solution is $$ x_p(t)=Φ(t)·\int_0^tΦ(s)^{-1}g(s)\,ds $$