  We propose the relaxation bootstrap method for the numerical solution of
multi-matrix models in the large $N$ limit, developing and improving the recent
proposal of H.Lin. It gives rigorous inequalities on the single trace moments
of the matrices up to a given "cutoff" order (length) of the moments. The
method combines usual loop equations on the moments and the positivity
constraint on the correlation matrix of the moments. We have a rigorous proof
of applicability of this method in the case of the one-matrix model where the
condition of positivity of the saddle point solution appears to be equivalent
to the presence of supports of the eigenvalue distribution only on the real
axis and only with positive weight. We demonstrate the numerical efficiency of
our method by solving the analytically "unsolvable" two-matrix model with
$\mathrm{tr}[A,B]^2$ interaction and quartic potentials, even for solutions
with spontaneously broken discrete symmetry. The region of values for computed
moments allowed by inequalities quickly shrinks with the increase of the
cutoff, allowing the precision of about 6 digits for generic values of
couplings in the case of $\mathbb{Z}_2$ symmetric solutions. Our numerical data
are checked against the known analytic results for particular values of
parameters.
