In this paper, we develop a first order (in time) numerical scheme for the binary fluid
surfactant phase field model. The free energy contains a double-well potential, a nonlinear coupling
entropy and a Flory-Huggins potential. The resulting coupled system consists of two Cahn-Hilliard
type equations. This system is solved numerically by finite difference spatial approximation, in
combination with convex splitting temporal discretization. We prove the proposed scheme is
unique solvable, positivity-preserving and unconditionally energy stable. In addition, an optimal
rate convergence analysis is provided for the proposed numerical scheme, which will be the first
such result for the binary fluid-surfactant system. Newton iteration is used to solve the discrete
system. Some numerical experiments are performed to validate the accuracy and energy stability
of the proposed scheme.