A computational framework of high order conservative finite difference methods to approximate the solution of a general system of N coupled nonlinear Schrödinger equations (N-CNLS) is proposed. Exact conservation of the discrete analogues of the mass and the system’s Hamiltonian is achieved by decomposing the original system into a sequence of smaller nonlinear problems, associated to each component of the complex field, and a modified Crank–Nicolson time marching scheme appropriately designed for systems. For a particular model problem, we formally prove that a method, based on the standard second order difference formula, converges with order τ+h2; and, using the theory of composition method, schemes of order τ2+h2 and τ4+h2 are derived. The methodology can be easily extended to other high order finite difference formulas and composition methods. Conservation and accuracy are numerically validated.