How to minimise two or more equations in R

I am trying to find an iterative way to solve these two M-estimater equations with two unknown parameters.

For each patient, we measure his blood pressure twice $Y_{i1}$ and $Y_{i2}$ and note his alcohol consumption $X_i$. We have given the following M-estimators and have proven these give unbiased results:

$$\sum\limits_{i=1}^n\sum\limits_{j=1}^2\big(Y_{ij}-\beta_0-\beta_1X_i\big)=0\quad\mbox{and } \sum\limits_{i=1}^n\sum\limits_{j=1}^2\big(Y_{i}-\beta_0-\beta_1X_i\big)X_i=0$$

(Where using OLS or maximum likelihood we assume all the measurements are independent).

I know it is possible to solve these analytically, but in case these two equations would be very complex, how do I solve these numerically in R?

Is there something like nlm for multiple equations?

Solution

What I should have done is minimise the norm of the two equations, and nlm is able to minimise a function with two unknown variables. This R script works fine:

y <- data_long$y
x <- data_long$x

f <- function(beta){
  temp_1 <- sum(y - beta[1] - beta[2] * x) 
  temp_2 <- sum(x*(y - beta[1] - beta[2] * x))
  sqrt(temp_1^2 + temp_2^2)
}
m_estimator <- nlm(f, c(0,0))

With data_long my simulated data, and beta my two estimates.