A two-flux model for light propagation in turbid media is proposed based on the differential and integral forms of Pomraning-Eddington approximation. The accuracy of the differential form depends on the choice of the weighting function and the boundary conditions, the thickness of the medium, and the scattering properties of the medium. In the integral form, the boundary conditions are imbedded, the range of validity of the approach is extended, and it is also able to deal with anisotropically scattering inhomogeneous media. Numerical results for partial fluxes and the albedos are calculated and compared with published calculations.