arXiv:2611.13419v1 Announce Type: new Abstract: We construct Neumann series of Bessel functions representations for solutions and their first derivatives of the equation $-y^{\prime\prime}(x)+(q(x)+2\rho p(x))y(x)=\rho^2y(x)$. The series arise from Legendre expansions of two transmutation operator kernels. The obtained series converge uniformly on horizontal strips in the complex plane of the parameter $\rho$. We develop a recursive integration procedure for computing the coefficients of these representations. The knowledge of the first coefficients allows one to recover the potentials $p$ and $q$ which is useful for solving inverse problems. For separated boundary conditions, we estimate the characteristic functions whose zeros determine the eigenvalues. Two numerical examples illustrate the practical applicability of the obtained series representations to the solution of spectral problems.