Need a numerical solution to simultaneous non-linear equations? The nleqslv package is just what you’re looking for! The coding required is minimal; just define the equations you want solved in a function, set some initial values, and let ‘er rip.
Here’s an example that uses the method of moments to estimate the parameters of a beta-binomial distribution.
The Home for Wayward Statisticians 