I have an exercise where I have to derive the both $w_i^{(m-1)}$ and $z_i^{(m-1)}$ from the iterative weighted least squared updating equation $b^{(m)} = \left( X^\top W^{(m-1)} X \right)^{-1} X^\top W^{(m-1)} z^{(m-1)}$ for the BeetleMortality data with the probit link $\phi$. Then I have to implement it in R. From what I understand, $w_i^{(m-1)}=\frac{1}{\text{Var}(Y)} (\frac{\partial \mu_i}{\partial \eta_i}) ^2 = \frac{\phi'(\eta_i)^{(m-1)}}{n\mu_i^{(m-1)}(1 - \mu_i^{(m-1)})}$ and $z_i^{(m-1)} = \eta_i^{(m-1)} (y_i -\mu_i) \frac{\partial \eta_i}{\partial \mu_i}= \eta_i^{(m-1)} + \frac{y_i + \mu_i^{(m-1)}}{\phi'(\eta_i)^{(m-1)}}$ where $\phi'$ is the normal distribution PDF. Now these could very well be wrong, but I like to think not. Now, my R implementation is as follows:

        
    library(glmx)
    library(msme)
    data("BeetleMortality")
    data <- as.data.frame(BeetleMortality)
    
    # initializing variables & coefficients
    x <- data$dose
    X <- cbind(1, x)  # design matrix with intercept
    y <- data$died
    n <- data$n
    b <- rep(0, ncol(X))
    print(X)
    
    # IRWLS algorithm
    for (iter in 1:5) {
      eta <- X %*% b
      mu <- pnorm(eta)  # fitted probabilities CDF

      # computing weights
      phi_eta <- dnorm(eta)  # ftandard normal PDF at eta    
      w <- (phi_eta^2) / (mu * (1 - mu))
      
      # computing working response
      z <- eta + (y - mu) / (phi_eta)
      
      # updating coefficients
      W <- diag(as.vector(w))
      b_new <- solve(t(X) %*% W %*% X) %*% (t(X) %*% W %*% z)
      b <- b_new  # update coefficients
      }
    
    cat("Estimated coefficients (IRWLS):\n")
    print(b)

    glm_fit <- glm(cbind(y, n - y) ~ dose, data = BeetleMortality, family = binomial(link = 
    probit))
    cat("Estimated coefficients (glm):\n")
    print(coef(glm_fit))
    
The code runs, but after only one iteration I get NaN values for the coefficients. I have tried re-deriving $w_i^{(m-1)}$ and $z_i^{(m-1)}$, but I always get the same. As you can see, I have a GLM at the end of the code, and runs perfectly fine. I am relatively new to all this so it could very well be that I made some glaring mistake which I cannot find. Any feedback is appreaciated!