General form of "Variation of Parameters" Recently I found out about the method to solve differential equations called "variation of parameters". However, I have only seen the form of this method for 2$^{nd}$ order DEs, namely
$$u_1'y_1+u_2'y_2=0$$
$$u_1'y_1'+u_2'y_2'=\frac {Q(x)}{a_2}$$
My question is:
Is there a general form for variation of parameters for the general n$^{th}$ order differential equation?
 A: A general form of variation of parameters works for first-order linear systems 
(and thus in particular $n$'th order linear DE's, which can be translated into first-order linear systems)
$$ X' = A(t) X + b(t)$$
Let $\Psi(t)$ be a fundamental matrix for the homogeneous system, i.e. 
solution of $\Psi' = A(t) \Psi$ with $\Psi(0) = I$.  Write 
$X(t) = \Psi(t) U(t)$.  Then the differential equation becomes
$$ \Psi' U + \Psi U' = A \Psi U + \Psi U' = A \Psi U + b $$
which simplifies to $\Psi U' = b$, so that $U = \int \Psi^{-1} b \; dt + const$.
A: Once you understand the basic idea of the method you can apply it to differential equations of any order. 

First Order Case
First consider the first order differential equation,
$$ a(t)\frac{d y}{dt}+b(t) y= F(t) ,$$
and assume that $y_{\mathrm h}$ solves the homogeneous case (i.e. $F(t)=0$).
We can write any function in the form of the product,
$$ \boxed{f(t) = A(t) \cdot y_{\mathrm h}(t)},$$
if we substitute this product into the differential equation we will be able to exploit the fact that $y_{\mathrm h}$ solves the homogeneous equation to make certain simplifications. 
$$a(t) \frac{d}{dt} \left( A(t) y_{\mathrm h}(t) \right) + b(t) \left( A(t) y_{\mathrm h}(t) \right) = F(t),$$
$$\Rightarrow a(t) \left( \frac{dA}{dt}  y_{\mathrm h}(t) + \color{blue}{A(t) \frac{d y_{\mathrm h}}{dt}}
  \right) + \color{blue}{b(t) A(t) y_{\mathrm h}(t) } = F(t),$$
$$\Rightarrow a(t) \frac{dA}{dt}  y_{\mathrm h}(t)  + A(t) \color{blue}{\left[ a(t)  \frac{d y_h}{dt}+b(t) y_{\mathrm h}(t) \right]} = F(t),$$
the blue expression is equal to $\color{blue}{0}$ by the definition of $y_{\mathrm{h}}$ leaving us with,
$$\boxed{a(t) \frac{dA}{dt}  y_{\mathrm h}(t)  = F(t)},$$
If we can solve this equation for $A$ then we will have found a particular solution, in the form of $f$, to the original differential equation. 

Example for 1st order case
As an example of the preceding section, consider the differential equation,
$$ \frac{dy}{dt}+y = \sin(t).$$
The homogeneous solution is $ y_{\mathrm{h}}(t) = y_0 \exp(-t)$. 
The differential equation for $A$ is then,
$$ \frac{dA}{dt} = e^t \sin(t),$$
so that we have,
$$\boxed{ A(t) = \int^t e^{t'} \sin(t') dt'} ,$$
and the particular solution to the original equation is given by,
$$ \boxed{f(t) = f_0 e^{-t} \int^t_0 e^{t'} \sin(t') dt'}. $$

Second Order Case
In the second order case we are solving,
$$ a(t) y'' + b(t) y' + c(t) y = F(t),$$
which generally has two linearly independent solutions ( $y^{(1)}_{\mathrm{h}}$ and $y^{(2)}_{\mathrm{h}}$).
We can then write any function as,
$$ f(t) = A(t) y^{(1)}_h + B(t) y^{(2)}_h(t),$$
substituting this in the differential equation will gives us differential equations for $A$ and $B$ just like in the first order case. For the third order case we write $f(t)$ as a linear combination of three linearly independent homogeneous solutions and the pattern continues for higher order differential equations. 
This answer isn't yet complete I'll be coming back later to fill in some details
A: Yes, there is a general form, valid also for linear equations whose coefficients are non-constant.
Method
Given an $n$-th order linear differential non homogeneous equation
$$
\sum_{k=0}^n a_k(x) y^{(k)}=b(x),\qquad a_n(x) \neq 0,
$$
and given $y_1(x),\ldots,y_n(x)$, $n$ independent solutions of the associated homogeneous equation, the following linear combination 
$$
y_p(x)=\sum_{h=1}^nc_h(x)y_h(x)
$$
is a solution if the coefficients $c_1(x),\ldots,c_n(x)$ satisfy the following linear system
$$
\left\{
\begin{align}
&\sum_{h=1}^nc'_h(x)y^{(k)}_h(x)=0, \qquad k=0,\ldots,n-2.\\
&\sum_{h=1}^nc'_h(x)y^{(n-1)}_h(x)=b(x),
\end{align}
\right.
$$
Solving the system and integrating the resulting $c'_1(x),\ldots,c'_n(x)$, provides the solution.
Explanation
Suppose you have a $n$-th order linear differential non homogeneous equation
$$
\sum_{k=0}^n a_k(x) y^{(k)}=b(x),\qquad a_n(x) \neq 0,
$$
and suppose you know $n$ independent solutions of the associated homogeneous equation $y_1(x),\ldots,y_n(x)$.
To find a solution of the non-homogeneous equation, try a solution like the following:
$$
y_p(x)=\sum_{h=1}^nc_h(x)y_h(x).\label{1}\tag{1}
$$
The first derivative is
$$
y'_p(x)=\sum_{h=1}^nc'_h(x)y_h(x)+\sum_{h=1}^nc_h(x)y'_h(x),
$$
but we suppose that the $c$'s are such that 
$$
\sum_{h=1}^nc'_h(x)y_h(x)=0
$$
i.e., they are non-constant, but behaves like constants in the whole, so that
$$
y'_p(x)=\sum_{h=1}^nc_h(x)y'_h(x).
$$
The second derivative is
$$
y'_p(x)=\sum_{h=1}^nc'_h(x)y'_h(x)+\sum_{h=1}^nc_h(x)y''_h(x),
$$
but, again, we suppose that the $c$'s are such that 
$$
\sum_{h=1}^nc'_h(x)y'_h(x)=0,
$$
so that
$$
y''_p(x)=\sum_{h=1}^nc_h(x)y''_h(x).
$$
Going on, we have in general
$$
y^{(k)}_p(x)=\sum_{h=1}^nc_h(x)y^{(k)}_h(x),\qquad k=0,\ldots,n-1,
$$
while, for $k=n$ we leave
$$
y^{(n)}_p(x)=\sum_{h=1}^nc'_h(x)y^{(n-1)}_h(x)+\sum_{h=1}^nc_h(x)y^{(n)}_h(x),
$$
with the additional conditions on the $c$'s
$$
\sum_{h=1}^nc'_h(x)y^{(k)}_h(x)=0, \qquad k=0,\ldots,n-2.\label{2}\tag{2}
$$
Substituting the expression for the derivatives in the differential equation, and remembering that the $y_h(x)$ are solutions of the homogeneous equation, all terms go away, it only remains
$$
\sum_{h=1}^nc'_h(x)y^{(n-1)}_h(x)=b(x),\label{3}\tag{3}
$$
so \eqref{1} is a solution if \eqref{2} and \eqref{3} are satisfied.
Solving this linear system with respect to $c'_h(x)$, and integrating, provides the solution of the non-homogeneous equation we looked for. 
