I need to calculate the estimators of a Probit Model, the thing is that I need to do it "manually", I can't use commands like lm or GLM. I already did it with those commands and I have the results, but now I need to do it with the formulas, maximize the Log-likelihood.

Hi,
You did the easy part. It's not clear to what level of external function one is allowed to use. You'll definitely need your matrix math (building a model matrix, error measurements, etc). You know, transposing, taking the inverse, matrix multiplication....solving a system of equations. Start by googling "design matrix", or "model matrix". To what level you're supposed to solve numerical integration or derivatives isn't clear from your post. No one would have the time to do that for you.

probit.nll <- function (beta) {
# linear predictor
eta <- W %*% beta
# probability
p <- pnorm(eta)
# negative log-likelihood
-sum((1 - Z) * log(1 - p) + Z * log(p))
}
probit.gr <- function (beta) {
# linear predictor
eta <- W %*% beta
# probability
p <- pnorm(eta)
# chain rule
u <- dnorm(eta) * (Z - p) / (p * (1 - p))
# gradient
-crossprod(W, u)
}
fit <- optim(beta, probit.nll, gr = probit.gr, method = "BFGS", hessian = TRUE)

I found this in another post somewhere, so if anyone needs to do a probit model manually this worked for me.