variation of parameters for systems
Once you can solve the unforced system x' = A x, variation of parameters is the universal method for adding any forcing g(t). Its idea is the same clever move that works for a single equation: take the homogeneous solution, which has a constant sitting in it, and let that constant become a function of time, then solve for what the function must be. You vary the parameter — hence the name.
Walk it through. The homogeneous solution is x_h = Phi(t) c, where Phi(t) is a fundamental matrix (for constant coefficients, e^(At)) and c is a constant vector. To solve x' = A x + g(t), guess that the answer has the same shape but with c promoted to a vector function: x = Phi(t) u(t). Differentiate by the product rule: x' = Phi' u + Phi u'. Since Phi' = A Phi, the first term is A Phi u = A x, so substituting into the equation, the A x parts cancel cleanly and what is left is Phi u' = g(t). That is the whole engine: solve for u' = Phi^(-1) g(t), then integrate, u(t) = integral of Phi^(-1)(s) g(s) ds. The particular solution is x_p = Phi(t) times that integral. The deliberate guess was engineered so that everything except the forcing cancels — just like an integrating factor is engineered to collapse a left side into a single derivative.
For constant-coefficient systems Phi(t) = e^(At) and Phi^(-1)(s) = e^(-As), so the formula becomes the matrix convolution x_p(t) = integral from t0 to t of e^(A(t - s)) g(s) ds. This is the Duhamel formula, and it is completely general: it asks nothing of g(t) beyond being integrable — no clever guessing, no special-form restriction. The price is that you must be able to integrate e^(-As) g(s), which may not have an elementary antiderivative; then the answer stands as an honest integral.
To solve x' = A x + g(t) with A = [0, 1; -1, 0] and g(t) = (0, cos t), use Phi(t) = e^(At) = [cos t, sin t; -sin t, cos t]. Form Phi^(-1)(s) g(s), integrate component by component, and multiply by Phi(t). The cancellation Phi u' = g happened automatically; you never had to guess the form of the answer, even though the forcing cos t is at the system's own frequency (where resonance makes a t-growth term appear in x_p).
Promote the constant vector c to a function u(t); the engineered cancellation leaves Phi u' = g, which you solve by one matrix inverse and one integral.
Variation of parameters is fully general but ends in an integral that may have no elementary form; undetermined coefficients is faster for nice forcing, yet only variation of parameters always works.