An imaginary-time nonuniform mesh method for diagonalizing multidimensional quantum Hamiltonians is proposed and used to find the first 50 eigenstates and energies of up to D=5 strongly interacting spinless quantum Lennard-Jones particles trapped in a one-dimensional harmonic potential. We show that the use of tailored grids allows exploitation of the symmetries of the system (in our case the D! degeneracy derived from all possible permutations of distinguishable particles), reducing drastically the computational effort needed to diagonalize the Hamiltonian. This leads to a favorable scaling with dimensionality, requiring for the five-dimensional system four orders of magnitude fewer grid points than the equivalent regular grid. Solutions to both bosonic and fermionic counterparts of this strongly interacting system are constructed, the bosonic case clustering as a Tonks-Girardeau crystal exhibiting the phenomenon of fermionization. The numerically exact excited states are used to describe the melting of this crystal at finite temperature.