In this case we will require deriving a new formula for variation of parameters for systems.  The derivation now will be much simpler than the when we first noticed variation of parameters.

First assume that X(t) be a matrix whose ith column is the ith linearly independent solution to the system,

x?' = A x?

Now this can be illustrated that X(t) will be a solution to the subsequent differential equation.

X' = AX    ......................(1)

It is nothing more than the original system along with the matrix instead of the original vector.

We are going to attempt and get a particular solution to,

x?' = A x? + g? (t)

We will suppose that we can get a solution of the form,

x?_p = X(t) v? (t)      

Here we will need to find out the vector v? (t). To do it we will require plugging this in the nonhomogeneous system. Keep in mind to product rule the particular solution whiler plugging the guess in the system.

X' v? + X v? = AX v? + g?

See that we dropped the "(t)" part of things to identify the notation a little. Here by using (1) we can rewrite this a little.

X' v? + X v? = X' v? + g?

X v? = g?

Since we formed X using linearly independent solutions we identify that det(X) should be nonzero and this in turn implies that we can get the inverse of X. Therefore, multiply both sides with the inverse of X.

v? = X-¹ g?

This time all that we require to do is integrate both sides to get v? (t).

v? (t) = ∫ (X-¹ g?) dt

When with the second order differential equation case we can ignore any constants of integration. So, the particular solution is,

x?_p = X ∫ (X^-1 v? (t)) dt