Summary: We study the thermodynamics of kinks in the presence of quasi-exactly-solvable (QES) potentials using a Langevin code implemented on a massively parallel computer. Specifically, we study the potential \(V(\varphi ) = \varphi ^{2} + (\lambda /\nu )/(1 + \nu \varphi ^{2})\) for which certain exact solutions of the Schrödinger equation are known provided the parameters \((\lambda \) and \(\nu )\) satisfy certain constraints. For \(\lambda 1\) this potential has two degenerate minima, otherwise it has only one minimum. Thus, a system with this potential is capable of describing a second-order phase transition (a la the \(\varphi ^{4}\) model; however, no exact solutions of the Schrödinger equation for this potential exist). We have obtained an exact solution for the kink in this model. In addition, we have calculated (in a 1D model) such quantities as the probability density function (PDF), field configuration and field-field correlation functions both above and below the transition temperature \((T_{c})\) for several temperatures. These quantities help us understand the contribution to the specific heat from coherent structures such as domain walls (kinks) as opposed to the contribution from lattice vibrations. We have calibrated our results against known exact solutions for limiting cases with very high accuracy.