In GLMs, why do we solve score(beta)=0 instead of just minimizing the negative log-likelihood? When we search for a numerical way to find $\hat{\beta}$ in a GLM (say, a logistic regression), we could do a numerical optimization (minimization) of the negative log-likelihood. 
But instead, we go one step further, and first compute the score function, to then equvialently look for the value of $\beta$ that sets $s(\beta) = 0$.
Why do we do this extra step? They are both iterative numerical procedures. Is solving this equation system computationally easier than minimizing the negative log-likelihood?
 A: The solutions to the score equation are critical points of the objective function in your optimisation, so generally the fitted coefficient estimator should solve the score equation.  This is not really a "step further" than using numerical techniques; it is just a way of characterising the fitted value, which is the point that numerical solutions should be moving towards.
Optimisation via IRLS: As a practical matter, numerical algorithms to fit GLMs generally use the default method of iteratively-reweighted least squares optimisation, which I think is equivalent to Fisher scoring for this class of models.  This is the default fitting-method in the R function glm (see comments on the method variable of the function in the documentation).  The estimation method can be derived via an argument involving linear approximation to the score equation, and linear approximation to the mean parameter in each iteration. 
The GLM is characteristed by having some exponential family distribution with scale parameter $\phi$ and mean parameter $\mu = h (\eta) = h(\boldsymbol{\text{x}} \beta)$ where the function $h$ is the canonical inverse-link and $\eta$ is a linear function of the coefficients of interest.  Using this distribution and taking a linear approximation to the score function yields:
$$s(\beta) = \frac{\partial \ell}{\partial \beta} (\beta) \approx \phi^{-1} \boldsymbol{\text{x}}^{\text{T}}(\boldsymbol{\text{y}} - \mu).$$
For some iteration value $\beta_{(k)}$ we can form the corresponding mean $\mu_{(k)} = h(\eta_{(k)}) = h(\boldsymbol{\text{x}} \beta_{(k)})$ and a linear approximation gives $\mu \approx \mu_{(k)} + \boldsymbol{\text{W}}_{(k)} \boldsymbol{\text{x}} (\beta - \beta_{(k)})$ with weight matrix $\boldsymbol{\text{W}}_{(k)} = \nabla_\eta h (\eta_{(k)})$.  Now, if we substitute this latter approximation into the approximation for the score equation we get:
$$\begin{equation} \begin{aligned}
s(\beta) = \frac{\partial \ell}{\partial \beta} (\beta) \approx \phi^{-1} \boldsymbol{\text{x}}^{\text{T}}(\boldsymbol{\text{y}} - \mu) 
&= \phi^{-1} \boldsymbol{\text{x}}^{\text{T}}(\boldsymbol{\text{y}} - \mu_{(k)} - (\mu - \mu_{(k)})) \\[8pt]
&\approx \phi^{-1} \boldsymbol{\text{x}}^{\text{T}}(\boldsymbol{\text{y}} - \mu_{(k)} - \boldsymbol{\text{W}}_{(k)} \boldsymbol{\text{x}} (\beta - \beta_{(k)})) \\[8pt]
&= \phi^{-1} \boldsymbol{\text{x}}^{\text{T}} \boldsymbol{\text{W}}_{(k)} ( \boldsymbol{\text{W}}_{(k)} ^{-1} (\boldsymbol{\text{y}} - \mu_{(k)}) + \boldsymbol{\text{x}} \beta_{(k)} - \boldsymbol{\text{x}} \beta ) \\[8pt]
&= \phi^{-1} \boldsymbol{\text{x}}^{\text{T}} \boldsymbol{\text{W}}_{(k)} ( \boldsymbol{\text{z}}_{(k)} - \boldsymbol{\text{x}} \beta ), \\[8pt]
\end{aligned} \end{equation}$$
where $\boldsymbol{\text{z}}_{(k)} = \boldsymbol{\text{W}}_{(k)} ^{-1} (\boldsymbol{\text{y}} - \mu_{(k)}) + \boldsymbol{\text{x}} \beta_{(k)}$ is the adjusted response, which is calculable from the $k$th iteration.  This approximation to the score function corresponds to a weighted linear regression.  Solving the score equation we obtain the next iteration for the coefficient vector, which has the form of a weighted-least-squares estimator:
$$\hat{\beta}_{(k+1)} = (\boldsymbol{\text{x}}^\text{T} \boldsymbol{\text{W}}_{(k)} \boldsymbol{\text{x}})^{-1} (\boldsymbol{\text{x}}^\text{T} \boldsymbol{\text{W}}_{(k)} \boldsymbol{\text{z}}_{(k)}).$$
This method is a simple variation on the method of least-squares estimation used in multiple linear regression with Gaussian errors.  Generally these algorithms start with a base iteration using least-squares estimation and then iterate the IRLS algorithm to within some tolerance of the true fitted value.  Note that even once the fitted value is obtained via IRLS, the score function is still useful, since it is related to the variance of the coefficient vector.
