The below are the estimating equations for the mean and dispersion parameters based on beta-binomial distribution.

Which, I captured from An empirical investigation of different operating characteristics of several estimators of the intraclass correlation in the analysis of binary data

I need to know, how can I solve two estimating equations with two unknown parameters in R?

(Note: This is not a GLM model)

I tried it with randomly generated beta-binomial variates using rootSolve R package, but didn't work...