Added a few lines so fitting correlation functions with no imaginary …

…part does not crash

qutip/solver/heom/bofin_baths.py

@@ -393,7 +393,7 @@ def correlation_function(self, t, **kwargs):
         ----------
             The correlation function as an array or float at time t
         """
-        if self.spectral_density(1) == -self.spectral_density(-1):
+        if np.isclose(self.spectral_density(t), -self.spectral_density(-t)):
 
             def integrand(w, t):
                 return (