It was already mentioned in a comment above about a topic in physics.SE, that it may be related with inverse scattering problem. I only may add, that there is precise method of construction of Shroedinger operator p2+V(q) with arbitrary finite spectrum for one-dimensional case. It is very well developed due to application to soliton theory. The polynomial operators V(q) have some problems due to infinities for q→±∞, but for V(q) with fast convergence to zero there are hundreds papers. Especially simple construction may be found, e.g., in M. Ablowitz, H. Segur, Solitons and the Inverse Scattering Transform, SIAM, Philadelphia, 1981.
The idea, that if you have some L1=p2+u(q) with given spectrum, there is an elegant method to add yet another eigenvalue ζ. You should consider equation L1g=ζg, find the solution g(q). Then L2=p2+u+2(ln(g))″ saves all eigenvalues of L1, but also has new eigenvalue ζ with eigenfunction 1/g(q). So you may start with u=0 and add eigenvalues one after another.
Yet, as I mentioned here is a problem with polynomial case, it is one-dimensional and may be only useful thing is idea to search for trick like g→1/g with adding eigenvalues...
[ADD] (1). I think, it is possible also to consider products Lk(qk), etc.
(2). It is not unique, e.g., isospectral deformation of L is described by KdV equation (yet, it was also mentioned in a reference in Physics.SE thread)
This post has been migrated from (A51.SE)