20個の2変数x,yデータで以下です。

ロジスティック回帰は一般化線形モデル（GLM）で扱うモデルのうち、従属変数が二項分布をとる場合についてロジット変換を使って回帰モデルを立てるもので、Rではglm()でfamilyをbinomial、リンク関数にlogitを使用します。
質問にあるモデル式ではリンク関数の記述が省略されているので、binomialの省略時デフォルトであるlogitが使用されます。

ans <- glm(y ~ x, d, family = binomial(logit))

GLMにおけるパラメータ推定は最尤法が用いられ、反復再重み付け最小二乗法（IRLS）が使用されます。最尤法の詳しい説明は参考書やGLMの解説をしたサイトなどいろいろあると思いますので、そういったものを参考にして下さい。
http://www7.atwiki.jp/hayatoiijima/sp/pages/18.html
乱暴に言いますと従属変数をロジット変換して線形回帰に持ち込むことで、データがロジスティック曲線にフィットするようなパラメータ推定を行えるようにしている、という感じでしょうか。

さてこのモデルにおいてどのような演算処理が行われているのかというご質問ですが、ロジスティック回帰の手順をRのコードで書いたものを付けておきますので参考にしてみて下さい。ここに書いたコードはこちらから拝借しました。
https://socserv.socsci.mcmaster.ca/jfox/Books/Companion-1E/Ch8-script.txt

まずデータフレームの読み込み

d <- data.frame(x = c(43.2808, 39.1153, 38.9163, 38.2939, 45.0221, 38.5792, 35.7308, 41.1779, 39.0925, 39.0529, 38.5518, 35.9808, 39.507, 40.0184, 39.4895, 42.835, 38.7532, 37.0281, 41.0135), y = c(1, 0, 0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 0, 0))

ロジスティック回帰（John Fox版）

lreg <- function(x, y, max.iter=10, tol=1E-6){
# x はモデル行列（説明変数のベクトル）
# y は反応ベクトル（0 または 1）
# max.iter は最大反復回数
# tol は収束基準
x <- cbind(1, x)
b <- b.last <- rep(0, ncol(x))
it <- 0
while (it <= max.iter){
p <- as.vector(1/(1 + exp(-x %*% b)))
V <- diag(p * (1 - p))
var.b <- solve(t(x) %*% V %*% x)
b <- b + var.b %*% t(x) %*% (y - p)
if (max(abs(b - b.last)/(abs(b.last) + 0.01*tol)) < tol) break
b.last <- b
it <- it + 1
}
if (it > max.iter) warning('maximum iterations exceeded')
list(coefficients=as.vector(b), var=var.b, iterations=it)
}
lreg(d$x, d$y)

上のコードの出力

$coefficients
[1] -82.219446   2.022379

$var
                     x
  1812.2569 -45.043403
x  -45.0434   1.120192

$iterations
[1] 7

glm()の出力と比較すると同じ係数と切片に収束していることがわかります。

ans <- glm(y ~ x, d, family = binomial(logit))
summary(ans)$coefficients

              Estimate Std. Error   z value   Pr(>|z|)
(Intercept) -82.219446  42.570413 -1.931375 0.05343665
x             2.022379   1.058386  1.910814 0.05602845