A closed form solution for the one-dimensional Schrödinger equation with a finite number of \(\delta \) -interactions \(\begin{aligned} {\textbf{L}}_{q,{\mathfrak {I}}_{N}}y:=-y^{\prime \prime }+\left( q(x)+\sum _{k=1}^{N}\alpha _{k}\delta (x-x_{k})\right) y=\lambda y,\quad 0<x<b,\;\lambda \in {\mathbb {C}} \end{aligned}\) is presented in terms of the solution of the unperturbed equation \(\begin{aligned} {\textbf{L}}_{q}y:=-y^{\prime \prime }+q(x)y=\lambda y,\quad 0<x<b,\;\lambda \in {\mathbb {C}} \end{aligned}\) and a corresponding transmutation (transformation) operator \({\textbf{T}}_{{\mathfrak {I}}_{N}}^{f}\) is obtained in the form of a Volterra integral operator. With the aid of the spectral parameter power series method, a practical construction of the image of the transmutation operator on a dense set is presented, and it is proved that the operator \({\textbf{T}}_{{\mathfrak {I}}_{N}}^{f}\) transmutes the second derivative into the Schrödinger operator \({\textbf{L}}_{q,{\mathfrak {I}}_{N}}\) on a Sobolev space \(H^{2}\) . A Fourier-Legendre series representation for the integral transmutation kernel is developed, from which a new representation for the solutions and their derivatives, in the form of a Neumann series of Bessel functions, is derived.
