A coupled potential finite volume method for simulation of three-dimensional marine controlled-source electromagnetic (CSEM) response in anisotropic formation is developed. To circumvent ill-conditioning and convergence problems, Maxwell's equations are reformulated into coupled scalar-vector potentials with Coulomb gauge and its complement by applying a Helmholtz decomposition to the electric field. Yee's staggered girds, finite volume averaging and interpolation techniques are used to make the Helmholtz equations discrete. The resulting sparse and complex linear system in large-scale models is solved by a direct solver PARDISO. In order to improve the accuracy of the near field results without significantly reducing the computational efficiency, a method using difference fields is proposed to reduce the source singularity effect of anisotropic formation. The anisotropic modeling examples show that marine CSEM response is predominantly sensitive to reservoir vertical resistivity, not to reservoir horizontal resistivity, provided that the reservoir are thin and high-resistive; but the marine CSEM response is sensitive to both horizontal and vertical resistivity of the overburden on top of the reservoir.