We consider the thermal softening of crystals due to anharmonicity. Self-consistent methods find a maximum temperature for a stable crystal, which gives an upper bound to the melting temperature. Previous workers have shown that the self-consistent harmonic approximation (SCHA) gives misleading results for the thermal stability of crystals. The reason is that the most important diagrams in the perturbation expansion around harmonic theory are not included in the SCHA. An alternative approach is to solve a self-consistent Dyson equation (SCA) for a selection of diagrams, the simplest being a (3+4)-SCA. However, this gives an unsatisfactory comparison to numerical results on the thermal and quantum melting of two-dimensional (2D) Coulomb-interacting particles (equivalent to vortex-lattice melting in two and three dimensions). We derive an improved self-consistent method, the two-vertex-SCHA, which gives much better agreement to the simulations. Our method allows for accurate calculation of the thermal softening of the shear modulus for 2D crystals and for the lattice of vortex-lines in type-II superconductors.