  The Heun functions have wide application in modern physics and are expected
to succeed the hypergeometrical functions in the physical problems of the 21st
century. The numerical work with those functions, however, is complicated and
requires filling the gaps in the theory of the Heun functions and also,
creating new algorithms able to work with them efficiently.
  We propose a new algorithm for solving a system of two nonlinear
transcendental equations with two complex variables based on the M\"uller
algorithm. The new algorithm is particularly useful in systems featuring the
Heun functions and for them, the new algorithm gives distinctly better results
than Newton's and Broyden's methods.
  As an example for its application in physics, the new algorithm was used to
find the quasi-normal modes (QNM) of Schwarzschild black hole described by the
Regge-Wheeler equation. The numerical results obtained by our method are
compared with the already published QNM frequencies and are found to coincide
to a great extent with them. Also discussed are the QNM of the Kerr black hole,
described by the Teukolsky Master equation.
