Variation of parameters and solutions of composite products of linear differential equations

Given a basis of solutions to k ordinary linear differential equations ℓ j [ y ] = 0 ( j = 1 , 2 , … , k ) , we show how Greenʹs functions can be used to construct a basis of solutions to the homogeneous differential equation ℓ [ y ] = 0 , where ℓ is the composite product ℓ = ℓ 1 ℓ 2 … ℓ k . The construction of these solutions is elementary and classical. In particular, we consider the special case when ℓ = ℓ 1 k . Remarkably, in this case, if { y 1 , y 2 , … , y n } is a basis of ℓ 1 [ y ] = 0 , then our method produces a basis of ℓ 1 k [ y ] = 0 for any k ∈ N . We illustrate our results with several classical differential equations and their special function solutions.