def stability(pyscf_mf):
    """
    Test wave function stability and re-optimize SCF.

    Args:
        pyscf_mf: PySCF mean field object (e.g. `scf.RHF()`)

    Returns:
        pyscf_mf: Updated PySCF mean field object
    """
    new_orbitals = pyscf_mf.stability()[0]
    new_1rdm = pyscf_mf.make_rdm1(new_orbitals, pyscf_mf.mo_occ)
    pyscf_mf = pyscf_mf.run(new_1rdm)

    return pyscf_mf