We propose a new algorithm for solving a system of two nonlinear transcendental equations with two complex variables based on the Muller algorithm. The two-dimensional Muller algorithm is tested on systems of different type and is found to work comparably to Newton's method and Broyden's method in many cases. The new algorithm is particularly useful in systems featuring the Heun functions whose complexity may make the already known algorithms not efficient enough or not working at all. In those specific cases, the new algorithm gives distinctly better results than the other two 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.