Using C₃₂ molecule as the skeleton of multipole moments, two types of substitutional co-doped fullerene derivatives C₂₈B₂N₂ and C₂₈B₂P₂ were designed with 16 isomers. The electronic properties, linear polarizability α and first hyperpolarizability β were studied by CAM-B3LYP method of density functional theory (DFT). The results show that the HOMO-LUMO gap of the doped molecules becomes smaller, and the values of α and β of C₂₈B₂P₂ are larger than those of C₂₈B₂P₂. Additionally, dipole molecules in these molecules have large β, while the octapolar molecules have smaller β. The structure with excellent second-order nonlinear optical (NLO) response characteristics is selected. The calculation results of time-dependent density functional theory (TD-DFT) method show that, compared with C₃₂, the response range of absorption spectra of all doped structures becomes wider, the maximum intensity becomes weaker, and the position of the maximum absorption wavelength may be red-shifted or blue-shifted. Based on the sum-over-states (SOS) method, the two-level or three-level formulas are used to explain the origin of large β in two co-doped structures, and verify the related electronic excitation type is π→π* excitation.