Without going too deeply into differential equation theory, one way to solve the equation is to diagonalize the matrix of coefficients, which can be done by a constant transformation that is invisible to the second derivative, once again solving differential equations in a single variable, setting out initial conditions, and transforming back to the coupled equations. Each eigenvector, which is called a normal mode, has its own time dependence, that of the overall system is their sum. An alternative would be to postulate exponential solutions, note that the second derivative multiplies the left hand vector by the square of the frequency, and treat discovering the frequencies as an eigenvalue problem. Either route requires diagonalizing the tridiagonal matrix of coefficients.
Using matrix algebra, and assuming that matrix calculus differs from scalar calculus mostly through the necessity to maintain the order of factors in a product and to write inverses where they belong instead of quotients, it is easy to describe the solution to forced motion once it is assumed that free motion is governed by the matrix exponential solution to a system of second order equations with constant coefficients.