  A spectral parameter power series (SPPS) representation for solutions of
Sturm-Liouville equations of the form
$$(pu')'+qu=u\sum_{k=1}^{N}\lambda^{k}r_{k}$$ is obtained. It allows one to
write a general solution of the equation as a power series in terms of the
spectral parameter $\lambda$. The coefficients of the series are given in terms
of recursive integrals involving a particular solution of the equation
$(pu_{0}')'+qu_{0}=0$. The convenient form of the solution provides an
efficient numerical method for solving corresponding initial value, boundary
value and spectral problems.
  A special case of the considered Sturm-Liouville equation arises in relation
with the Zakharov-Shabat system. We derive an SPPS representation for its
general solution and consider other applications as the one-dimensional Dirac
system and the equation describing a damped string. Several numerical examples
illustrate the efficiency and the accuracy of the numerical method based on the
SPPS representations which besides its natural advantages like the simplicity
in implementation and accuracy is applicable to the problems admitting complex
coefficients, spectral parameter dependent boundary conditions and complex
spectrum.
