Abstract: We developed a high performance algebraic solver for nonlinear systems discretized from two-dimensional energy equations with three temperatures by a nine point scheme.The main idea is to solve the system by an inexact Newton method and preconditioned Krylov subspace methods in the frame of PNK and JFNK methods.Numerical experiments show the efficiency of the algebraic solvers.It is shown that our PNK method is 6 times faster than the nonlinear block Gauss-Seidel method. The JFNK and PNK methods are also compared.