  We give explicit formulas for a pair of linearly independent solutions of
$(py')'(x)+q(x)=(\lambda_1r_1(x)+\cdots+\lambda_dr_d(x))y(x)$, thus
generalizing to arbitrary $d$ previously known formulas for $d=1$. These are
power series in the spectral parameters $\lambda_1,\dots,\lambda_d$ (real or
complex), with coefficients which are functions on the interval of definition
of the differential equation. The coefficients are obtained recursively using
indefinite integrals involving the coefficients of lower degree. Examples are
provided in which these formulas are used to solve numerically some boundary
value problems for $d=2$, as well as an application to transmission and
reflectance in optics.
