Ab initio spin-flip configuration interaction (SF-CI) methods with the finite-field (FF) scheme are applied to the calculation of static second hyperpolarizabilities (γ) of several singlet diradical systems, i.e., the model H2 molecule under dissociation, p-quinodimethane, o-quinoid five-membered ring, and 1,4-bis(imidazole-2-ylidene)cyclohexa-2,5-diene (BI2Y) models. The SF-CI method using the UHF reference wave function provides the qualitatively correct diradical character (γ) dependence of γ in a wide range of a diradical character region for H2 under dissociation ahd p-quinodimethane as well as o-quinoid five-membered ring models. For BI2Y, which is a real diradical system, a non-negligible spin contamination is found in the spin-unrestricted Hartree-Fock (UHF) triplet state, which results in overestimations (SF-CIS) or underestimations (SF-CIS(D)) of γ. Such deficiencies are significantly reduced when using the pure spin state, i.e., the restricted open-shell HF (ROHF) triplet wave function as the reference wave function. These results indicate the applicability of the FF-SF-CI method starting with a pure or a nearly pure high-spin state to provide qualitative or semiquantitative γ for large-size diradical systems. For selected systems, these SF-CI results are also compared to the SF equation of motion coupled cluster singles and doubles (SF-EOM-CCSD) and to SF time-dependent density functional theory (SF-TDDFT) schemes. In particular, large amounts of Hartree-Fock exchange in the functional are required to obtain qualitatively correct dependence of γ on y in the case of p-quinodimethane.