The exact solution of spherical mean-field plus multi-pair interaction model with two non-degenerate j-orbits, which is an extension of the widely used standard (two-body) pairing model, is derived based on the Bethe–Richardson–Gaudin approach. The Bethe–Richardson–Gaudin equations in determining eigenstates and the corresponding eigen-energies of the model are provided and exemplified with up to three-pair interactions. With a suitable parameterization of the overall multi-pair interaction strengths, the model with one adjustable parameter and valence nucleons confined in the $$1d_{5/2}$$ and $$0g_{7/2}$$ orbits is applied to fit binding energies of $$^{102{\text {--}}112}$$Sn. It is shown that the ground-state occupation probabilities of nucleon pairs calculated from this model and those from the standard pairing model are almost the same with perfect ground-state overlap of the two models. A noticeable feature of the multi-pair interactions is that the even–odd staggering of pairing interaction strength appearing in the standard pairing model due to the Pauli-blocking is suppressed. As the result, the pairing interaction strength of the model only depends on the number of valence nucleon pairs in the system.