In this letter, a new algorithm for joint diagonalization of a set of matrices by congruence is proposed to compute the nonnegative joint diagonalizer. The nonnegativity constraint is imposed by means of a square change of variables. Then we formulate the high-dimensional optimization problem into several sequential polynomial subproblems using LU matrix factorization. Numerical experiments on simulated matrices emphasize the advantages of the proposed method, especially in the case of degeneracies such as for low SNR values and a small number of matrices. An illustration of blind separation of nuclear magnetic resonance spectroscopy confirms the validity and improvement of the proposed method.