ey = resid(lm(y ~ x2 + x3)) calculates the variance of y not explained by the regressors $x_2$ and $x_3$. On the other hand, ex = resid(lm(x ~ x2 + x3)) gives you the variance of $x$ not explained by $x_2$ and $x_3$.
Therefore, regressing ey over ex will calculate the
$\text{variance of }y \text{ not explained by }x_2\text{ and} x_3$
explained by the
$\text{variance of }x_1 \text{ not explained by }x_2\text{ and} x_3$
So you have eliminated the contribution of $x_2$ and $x_3$, and you are really just calculating the variance of $y$ explained by $x_1$, after having regressed both variables (dependent an independent) over $x_2$ and $x_2$ to effectively eliminate the effect of these two additional explanatory variables.
Now you just have to keep in mind that the regression through the origin with one single regressor is $\hat \beta=\frac{\displaystyle\sum_{i=1}^{n}Y_iX_i}{\displaystyle\sum_{i=1}^{n}X_i^2}$, and that the $-1$ in coef(lm(ey ~ ex - 1)) is the call for the OLS without intercept.
The reason why the intercept is eliminated is because the intercept is a regressor in its own right. Even if we didn't specify it in the call for ey and ex the regression over $1$ was there in the model matrix.
The package {swirl} in R contains step-by-step practice in this technique of picking one regressor, replacing all other variables by the residuals of their regressions against that one.