It is a standard approach to consider the maximum likelihood estimation procedure for the estimation of parameters in statistical modelling. The sample likelihood function has a closed form representation only if the two densities in the integrand are conjugate to each other. In case of any non- conjugate pair, no closed form representation exists. In such situations, we need to approximate the integral by making use of some numerical techniques. A first or second order Laplace approximation or the (adaptive) Gauss-Hermite quadrature method can be applied in order to get an approximative objective function. The resulting approximation of the likelihood function still needs to be numerically maximized with respect to all unknown parameters. For such a numerical maximization, all required derivatives are provided in the scope of this work. We explore the use of the (adaptive) Gauss-Hermite quadrature for Generalized Linear Mixed Models, when the conditional density of the response given the random effects is a member of the linear exponential family and the random effects are Gaussian.