In this article, the Galerkin method, based on quintic B-spline function as the shape and weight functions is described for the numerical solution of the second order coupled nonlinear Schrödinger equations. Finite difference and Crank–Nicolson schemes are used to discretize the time derivative and nodal parameters respectively. Three numerical problems are presented to assess the accuracy and capability of the proposed method. The maximum errors, norms and conserved quantities are calculated. The obtained numerical results show that the present scheme with higher order B-spline as basis and weight functions performs well and accurately. The numerical results are compared with analytical and published results.