Algorithm for estimating function maximizing parameters
You have a function $f:\mathbb{R}^3 \to \mathbb{R}$, which accepts three real numbers as input and produces a real number as output.  You want to find values $\omega, \alpha, \beta$ that maximize $f(\omega, \alpha, \beta)$.

This is a standard problem in mathematical optimization, and there are a variety of techniques for solving this problem.


Black-box methods. The methods that are easiest to use treat $f$ as a complete black box: all we need is  a way to evaluate $f$ at any input of our choice.  Example methods include hill climbing, coordinate descent, (well, in this case it will be coordinate ascent, since you want to maximize $f$), gradient descent (well, in this case it will be gradient ascent, since you want to maximize $f$), Newton's method, conjugate gradient methods, and more.
Black-box, plus derivatives. Some of these methods can be sped up significantly if we have an efficient way to evaluate each of the partial derivatives of $f$ at any input of our choice.  Sometimes, if $f$ has a simple representation, we can compute the partial derivatives of $f$ algebraically, which then gives us a way to evaluate these functions at any desired input.  The reason this speed things up is that some of the black-box methods effectively end up approximating the partial derivatives or 2nd partial derivatives (e.g., the gradient or Hessian matrix) as part of the process; if you have a faster way to compute them directly, then those methods become more efficient and more effective.
Dedicated methods. There are other methods available if the function $f$ has a particular algebraic form, e.g., linear programming, convex optimization, and more.  These probably won't apply to you, so they can be safely ignored in your situation.


Many tools, like MathLib, include built-in libraries to apply these methods for you.

In your case, the function $f$ is given by $f(\omega, \alpha, \beta) = \sum_{i=1}^m{\left\{ -\ln(v_i) - \frac{r_i^2}{v_i} \right\}}$ where $v_i = \omega + \alpha r_{i-1}^2 + \beta v_{i-1}^2$.  But those details aren't important; what's important is that, given $\omega, \alpha, \beta$, you can evaluate $f$ to learn the value of $f(\omega, \alpha, \beta)$.  That allows you to apply black-box optimization methods.

One generic issue to be aware of is that many of these methods may only produce a local optimum, not a global optimum.  There are a variety of heuristics for dealing with this, but generally none that are guaranteed to work in all cases.  Often, the methods are fairly sensitive to the choice of starting point, so if you have a good guess at a value of $\omega, \alpha, \beta$ that is likely to be near the optimal value, then starting the methods from there can be helpful.  Also, it can be helpful to apply the methods multiple times, from different (random) starting points and take the best observed over all iterations.