Express how additional covariates in a multiple regression may affect the coefficient of a variable in a simple linear regression In a simple linear regression
$$y=\alpha+\beta x + \varepsilon,$$
if $x$ is binary we can show that the slope coefficient ($\beta$) is equal to the difference in the outcome variable between levels of the predictor, i.e.
$$E(y\mid x) = \alpha+\beta x$$ so
$$E(y\mid x=0) = \alpha$$ and $$E(y\mid x=1) = \alpha+\beta$$.
We then have:
$$E(y\mid x=1) - E(y\mid x=0) = \beta.$$
However, in a multiple regression this usually is not the case. When we control other variables, some of the variation in $y$ originally explained by $x$ may be partitioned into the new covariate, which changes the estimate on $\beta$ so that the original relationship of $E(y\mid x=1) - E(y\mid x=0) = \beta$ will no longer be true.
I understand this intuitively and can show it with data, but I'm not sure how to express it mathematically with formulas (?).
An illustration using data in R:
set.seed(1)
# generate data
outcome    <- rnorm(n=100, mean = 10, sd = 1)
original_x <- rbinom(n=100, 1, prob=.5)
new_x      <- rnorm(n=100, mean = 10, sd = 10)

# Simple Linear Regression
coef(glm(outcome~original_x))[2]  # 0.08680985
mean(outcome[original_x==1])-mean(outcome[original_x==0]) # 0.08680985

# Multiple Regression
coef(glm(outcome~original_x + new_x))[2] # 0.1112377 

 A: Going to change up the notation a bit here, to match what you really did in your code.  We are assuming that $E(Y|X=x) = x\beta$, so, given $X = x$,
$$Y = x\beta + \epsilon$$
where $\beta = (\alpha, \beta_1)^T$, and $x$ is an $n \times 2$ matrix with observations as rows. Implicitly, we are also assuming that $x$ is full-rank.
As I'm sure you know, the OLS estimate for $\beta$ is
$$\hat{\beta}_{OLS} = (x^Tx)^{-1}x^Ty$$
where $y$ is our observed response. What happens if we add a new regressor? Let the vector of observations for this new regressor be $x_a$. Then, we can define a new design matrix $x_{new} = \begin{pmatrix} x & x_a\end{pmatrix}$ which is $n \times 3$. Our new regression coefficients are
$$\hat{\beta}_{new} = (x_{new}^Tx_{new})^{-1}x_{new}^Ty = \begin{pmatrix} x^Tx & x^Tx_a\\ x_a^Tx & x_a^Tx_a\end{pmatrix}^{-1}\begin{pmatrix}x^Ty \\ x_a^Ty \end{pmatrix}$$
where we're now estimating three coefficients, namely $\beta = (\alpha, \beta_1, \beta_2)$. Notice, that the cross terms $x^Tx_a$ and $x_a^Tx$ are in a sense "responsible" for why the estimates for $\alpha$ and $\beta_1$ change from $\hat{\beta}$ to $\hat{\beta}_{new}$.
From here, I'd like to make a qualification to your statement. In most cases, you are right; your coefficient estimates change. In the case when the columns of $x$ are orthogonal to $x_a$, however, the vectors $x^Tx_a$ and $x_a^Tx$ are both matrices of different sizes, with all elements zero. Further algebraic trickery (taking the inverse of a $2 \times 2$ block matrix, which I won't do here) will reveal that
$$\hat{\beta}_{new} = \begin{pmatrix} (x^Tx)^{-1}(x^Ty)\\  (x_a^Tx_a)^{-1}x_a^Ty \end{pmatrix}$$
Therefore, in this case, each regressor "stays in its own lane"; thus, the addition or removal of a certain regressor will not change the estimates for the others.
