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.