Mathematical and numerical analysis has been undertaken for a pair of coupled Cahn-Hilliard equations with a logarithmic potential and with homogeneous Neumann boundary conditions. This pair of coupled equations arises in a phase separation model of thin film of binary liquid mixture. Global existence and uniqueness of a weak solution to the problem is proved using Faedo-Galerkin method. Higher regularity results of the weak solution are established under further regular requirements on the initial data. Further, continuous dependence on the initial data is presented. Numerically, semi-discrete and fully-discrete piecewise linear finite element approximations to the continuous problem are proposed for which existence, uniqueness and various stability estimates of the approximate solutions are proved. Semi-discrete and fully-discrete error bounds are derived where the time discretisation error is optimal. An iterative method for solving the resulting nonlinear algebraic system is introduced and linear stability analysis in one space dimension is studied. Finally, numerical experiments illustrating some of the theoretical results are performed in one and two space dimensions.