diff --git a/GPy/core/gp.py b/GPy/core/gp.py
index 6441561b..2cff4341 100644
--- a/GPy/core/gp.py
+++ b/GPy/core/gp.py
@@ -70,7 +70,7 @@ class GP(Model):
         self.posterior, self._log_marginal_likelihood, grad_dict = self.inference_method.inference(self.kern, self.X, self.likelihood, self.Y, Y_metadata=self.Y_metadata)
         self.likelihood.update_gradients(np.diag(grad_dict['dL_dK']))
         self.kern.update_gradients_full(grad_dict['dL_dK'], self.X)
-        
+
     def log_likelihood(self):
         return self._log_marginal_likelihood
 
@@ -186,7 +186,7 @@ class GP(Model):
         """
         assert "matplotlib" in sys.modules, "matplotlib package has not been imported."
         from ..plotting.matplot_dep import models_plots
-        models_plots.plot_fit_f(self,*args,**kwargs)
+        return models_plots.plot_fit_f(self,*args,**kwargs)
 
     def plot(self, *args, **kwargs):
         """
@@ -207,7 +207,7 @@ class GP(Model):
         """
         assert "matplotlib" in sys.modules, "matplotlib package has not been imported."
         from ..plotting.matplot_dep import models_plots
-        models_plots.plot_fit(self,*args,**kwargs)
+        return models_plots.plot_fit(self,*args,**kwargs)
 
     def _getstate(self):
         """
diff --git a/GPy/core/model.py b/GPy/core/model.py
index 6a6fe1ba..710c1b22 100644
--- a/GPy/core/model.py
+++ b/GPy/core/model.py
@@ -339,7 +339,7 @@ class Model(Parameterized):
                     print "No free parameters to check"
                     return
 
-            gradient = self.objective_function_gradients(x)
+            gradient = self.objective_function_gradients(x).copy()
             np.where(gradient == 0, 1e-312, gradient)
             ret = True
             for nind, xind in itertools.izip(param_index, transformed_index):
@@ -350,7 +350,12 @@ class Model(Parameterized):
                 f2 = self.objective_function(xx)
                 numerical_gradient = (f1 - f2) / (2 * step)
                 if _debug:
+                    for p in self.kern.flattened_parameters:
+                        p._parent_._debug=True
                     self.gradient[xind] = numerical_gradient
+                    self._set_params_transformed(x)
+                    for p in self.kern.flattened_parameters:
+                        p._parent_._debug=False
                 if np.all(gradient[xind]==0): ratio = (f1-f2) == gradient[xind] 
                 else: ratio = (f1 - f2) / (2 * step * gradient[xind])
                 difference = np.abs((f1 - f2) / 2 / step - gradient[xind])