We present WarmSPy, a numerical code in Python designed to solve for the perturbations' equations in warm inflation models and compute the corresponding scalar power spectrum at CMB horizon crossing. In models of warm inflation, a radiation bath of temperature T during inflation induces a dissipation (friction) rate of strength Q ∝ T^c/ϕ^m in the equation of motion for the inflaton field ϕ. While for a temperature-independent dissipation rate (c = 0) an analytic expression for the scalar power spectrum exists, in the case of a non-zero value for c the set of equations can only be solved numerically. For c > 0 (c < 0), the coupling between the perturbations in the inflaton field and radiation induces a growing (decaying) mode in the scalar perturbations, generally parameterized by a multiplicative function G(Q) which we refer to as the scalar dissipation function. Using WarmSPy, we provide an analytic fit for G(Q) for the cases of c = {3,1,-1}, corresponding to three cases that have been realized in physical models. Compared to previous literature results, our fits are more robust and valid over a broader range of dissipation strengths Q ∈ [10^-7,10⁴]. Additionally, for the first time, we numerically assess the stability of the scalar dissipation function against various model parameters, inflationary histories as well as the effects of metric perturbations. As a whole, the results do not depend appreciably on most of the parameters in the analysis, except for the dissipation index c, providing evidence for the universal behaviour of the scalar dissipation function G(Q).