Abstract We apply the iteration of source function (IOSF) philosophy to the successive order of scattering method for solving the vector radiative transfer equation in the coupled atmosphere and ocean system. A major class of radiative transfer solvers only provides the radiation field at discrete viewing zenith angles. The radiation field at other angles is found by interpolation. The iteration of source matrix method integrates the product of the radiation field and source matrix at quadrature points to obtain the radiation field at arbitrary viewing angles. The resultant solution includes the radiation contributions from all scattering orders higher than one. The analytical single scattering solution is then added to find the total radiation field. The proposed scheme includes the benefits of both the IOSF interpolation and the analytical single scattering solution. Boundary conditions for a flat air–sea interface are fully considered. A test case of a coupled atmosphere and ocean system shows that this combined method improves the polarized radiation field greatly in comparison with the regular polynomial interpolation method.