So an estimator for for $\lambda$ using the first method of moment is simply the sample mean:
$$\mu_1 = E(X) = \lambda = \bar{X} = \hat{\mu_2}$$
For the second moment,
$$\operatorname{Var}(X) = E(X^2) - E(X)^2$$ $$\lambda = E(X^2) - \lambda^2$$ $$\mu_2 = E(X^2) = \lambda + \lambda^2$$ $$s^2 = \frac{1}{n}\sum_{i=1}^n (X_i - \bar{X})^2$$ $$s^2 = \frac{1}{n}\sum_{i=1}^n X_i^2 - \bar{X}^2$$ $$\frac{1}{n}\sum_{i=1}^n X_i^2 = s^2 + \bar{X}^2$$ $$\hat{\mu_2} = s^2 + \bar{X}^2$$ Making $\mu_2 = \hat{\mu_2}$, we get:
$$\mu_2 = \lambda + \lambda^2 = s^2 + \bar{X}^2 = \hat{\mu_2}$$ $$ \lambda + \lambda^2 = s^2 + \bar{X}^2$$
From here I've no idea how to factorise it so I have $\lambda$ all on one side.
So you have three potential estimators:
in the R simulation below, SDs of estimators are given instead of their variances. With $m = 100\,000$ iterations we can expect about two significant digits of accuracy.