After thinking it over, this discrepancy makes perfect sense.
Because CCSD (and MP2) use Hartree-Fock orbitals, we can easily separate the Lagrangian into a Hartree-Fock energy functional and a correlation energy functional. To get the orbital unrelaxed dipole moment, you want to compute the derivative of the Lagrangian energy with respect to an external electric field, except removing the orbital relaxation term in the correlation energy functional. You tried to accomplish this by eliminating all dependence of the Hartree-Fock orbitals and energies on the electric field. Unfortunately, this eliminates the dependence of the Hartree-Fock energy functional on your electric field, so you have completely eliminated the Hartree-Fock electronic dipole contribution. Unlike orbital relaxation, that is a very large contribution and goes far beyond the orbital unrelaxed approximation. As evidence, if you look at the breakdown of the final energy at the end of your unrelaxed output files, you see the final HF energy matches the final HF energy you computed with no perturbation.
You are not numerically differentiating the quantity you think you are. I can help with using Psi4, but devising the equations to perform this numerical differentiation for you is past what I have time to do.