We suggest a new computer-assisted approach to the development of turbulence theory. It allows one to impose lower and upper bounds on correlation functions using sum-of-squares polynomials. We demonstrate it on the minimal cascade model of two resonantly interacting modes when one is pumped and the other dissipates. We show how to present correlation functions of interest as part of a sum-of-squares polynomial using the stationarity of the statistics. That allows us to find how the moments of the mode amplitudes depend on the degree of nonequilibrium (analog of the Reynolds number), which reveals some properties of marginal statistical distributions. By combining scaling dependence with the results of direct numerical simulations, we obtain the probability densities of both modes in a highly intermittent inverse cascade. As the Reynolds number tends to infinity, we show that the relative phase between modes tends to $\pi/2$ and $-\pi/2$ in the direct and inverse cascades, respectively, and derive bounds on the phase variance. Our approach combines computer-aided analytical proofs with a numerical algorithm applied to high-degree polynomials.