An efficient method is proposed for evaluating first- and second-order nonadiabatic matrix elements of the form i(q;Q)j(q;Q)/Q and i(q;Q)2j(q;Q)/Q2, where i(q;Q) and j(q;Q) denote multiconfigurational self-consistent-field electron wave functions. The method is based on a finite-difference procedure and requires the numerical computation of symmetric overlaps of the type i(q,Q0-x)j(q,Q0 +x). It gives an accuracy which is quadratic in the nuclear displacement x for both the first- and second-order nonadiabatic coupling constants. The wave functions are separately optimized for each state and obtained through the direct second-order MCSCF method. The biorthogonal scheme of Malmquist is implemented that expresses i(q;Q) and j(q:Q) in an orthogonal common basis. The method is applied for the calculation of the nonadiabatic coupling elements and the Born-Oppenheimer corrections to the two lowest +2 states of NaLi+, relevant for analyzing the asymmetric charge exchange in the ion-atom collision Na+Li+Na++Li.. AE.