Schrödinger equation with finitely many $\delta$-interactions: closed form, integral and series representations for solutions

A closed form solution for the one-dimensional Schrödinger equation with a finite number of $\delta$-interactions $${\mathbf{L}}_{q,{\mathfrak{I}}_N}y:=-y^{\mathrm{\prime }\mathrm{\prime }}+(q(x)+\sum _{k=1}^N{\alpha }_k\delta (x-x_k))y=\lambda y,0<x<b,\lambda \in \mathbb{C}$$ is presented in terms of the solution of the unperturbed equation $${\mathbf{L}}_{q}y:=-y^{\mathrm{\prime }\mathrm{\prime }}+q(x)y=\lambda y,0<x<b,\lambda \in \mathbb{C}$$ and a corresponding transmutation operator ${\mathbf{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 ${\mathbf{T}}_{{\mathfrak{I}}_N}^f$ transmutes the second derivative into the Schrödinger operator ${\mathbf{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.