A new numerical scheme for the solution of liquid state integral equations using the Baxter factorization of the Ornstein-Zernike equation is proposed. For short range potentials the method yields reliable results over the whole fluid region, including the vicinity of the critical point, and opens up new possibilities for numerical study of the critical behavior of integral equation approximations. To demonstrate the effectiveness of the method, numerical results are compared with the analytical solution of the mean spherical approximation for a hard-core plus Yukawa tail interaction potential. Accurate results for the critical exponents δ, γ, and η for this model are obtained.