For many nonlinear physical systems, approximate solutions are pursued by conventional perturbation theory in powers of the non-linear terms. Unfortunately, this often produces divergent asymptotic series, collectively dismissed by Abel as "an invention of the devil." An alternative method, the self-consistent expansion (SCE), has been introduced by Schwartz and Edwards. Its basic idea is a rescaling of the zeroth-order system around which the solution is expanded, to achieve optimal results. While low-order SCEs have been remarkably successful in describing the dynamics of non-equilibrium many-body systems (e.g., the Kardar-Parisi-Zhang equation), its convergence properties have not been elucidated before. To address this issue we apply this technique to the canonical partition function of the classical harmonic oscillator with a quartic $gx^{4}$ anharmonicity, for which perturbation theory's divergence is well-known. We obtain the $N$th order SCE for the partition function, which is rigorously found to converge exponentially fast in $N$, and uniformly in $g\ge0$. We use our results to elucidate the relation between the SCE and the class of approaches based on the so-called "order-dependent mapping." Moreover, we put the SCE to test against other methods that improve upon perturbation theory (Borel resummation, hyperasymptotics, Pad\'e approximants, and the Lanczos $\tau$-method), and find that it compares favorably with all of them for small $g$ and dominates over them for large $g$. The SCE is shown to successfully capture the correct partition function for the double-well potential case, where no perturbative expansion exists. Our treatment is generalized to the case of many oscillators, as well as to any nonlinearity of the form $g|x|^{q}$ with $q\ge0$ and complex $g$. These results allow us to treat the Airy function, and to see the fingerprints of Stokes lines in the SCE.