How to obtain the eigensystem for the quantum dipole chain? Lets consider a quantum dipole Hamiltonian
\begin{equation}
H=\sum_{i=1}^{L}\frac{p_i^2}{2m_i} + \frac{1}{2}k_i x_i^2 + \frac{\Lambda}{2}\sum_{i,j=\langle n,n\rangle}\frac{1}{|x_i - x_j|^3},
\end{equation}
where $\langle n,n\rangle$ stands for nearest neighbors, i.e., the dipole interaction is taken over the nearest neighbor pairs only, and the Hamiltonian also includes a harmonic potential.
I don't think there is an analytic way to obtain the eigenvalues and eigenvectors of the above Hamiltonain. Hence, what would be the best numerical way? 
A simple way is to use the harmonic oscillator energy basis that will result in a non-diagonal quantum dipole Hamiltonian matrix and diagonalize it using standard routines. The problem with this approach is that even to obtain the first few energy levels and the corresponding eigenvectors accurately the truncation of basis for each harmonic oscillator needs to be very large ~$20$. Hence going to large $L$ is problematic due to the sharp increase in Hilbert space dimension, only allowing for a maximum of $L=4$.
Another approach with similar problems could possibly be to represent the Hamiltonain in the position basis.
Or is there some other alternative I haven't looked at that would simplify this problem and help go beyond $L=4$?
 A: In short: not really, no. Many-body quantum mechanical problems are just Hard, and if what you want is the full solution then the problem is exponentially hard in $L$; you might be able to delay the onset of the difficulties by a bit by being clever about things, but the gains in $L$ will be limited.
The "being clever about things" is normally a matter of choosing a basis which is suited well to the problem at hand, and for your configuration the harmonic-oscillator basis is unlikely to be a one of those, because the basis takes no account of the actual hamiltonian of the problem. To do that, one thing you could do is to (i) solve the classical problem to get the classical equilibrium positions, and then (ii) choose a set of displaced number states at the classical equilibrium positions (which then (iii) need to be orthogonalized, ideally in a symmetric fashion), and use those as a basis. Or a bunch of other possible schemes to use the particular features of the problem to cut down on the basis size.
Ultimately, your problem isn't really that different from that of having $L$ different electrons in an atom, subject to a central potential and interacting with each other, and that is also, in general, an exponentially hard problem. How do we deal with this? By being very drastic in the (irrelevant) parts of the Hilbert space that we cut out of the calculation, in a two-stage approach:


*

*First, by using a self-consistent-field approach to solve for the optimal single-particle basis, normally known as the Hartree-Fock method. Here the idea is that you solve for the ground state by pretending that all the electrons are the same one, so you start with some decent guess for the ground state, use that to calculate the mean-field potential, solve for the ground state, put that back in as a mean-field potential, and iterate until it (hopefully) converges, possibly followed by a variational approach where you vary any relevant parameters until the energy is as low as possible.
This will get you a Hartree-Fock basis, i.e. a set of single-particle orbitals that incorporate as much as possible of the diagonal, uncorrelated part of the inter-particle interaction.

*After that, you take a bunch of multi-particle configurations using the HF single-particle orbitals. You can just take all of them and have something that is formally exact; this will still have an exponential dimension, but the improvement in the basis might be enough to push to a reasonable number of particles.
If you want to get serious about this, though, you normally use one of multiple post-Hartree-Fock methods, which essentially try to cut out configurations which are not essential to the description, while extending the basis in the directions that do the most good. This class of methods is capable of capturing some quite correlated eigenstates for modest computational resources, but they do require a lot of coding to get right.
Generally speaking, though, the many-body problem is complicated, and though each branch of QM where it appears has its ways of mitigating the ensuing mess, they are all just fighting a polynomial fight against an exponential scaling.
