diff --git a/GPy/models/GP.py b/GPy/models/GP.py
index c640e529..2afa4252 100644
--- a/GPy/models/GP.py
+++ b/GPy/models/GP.py
@@ -133,7 +133,7 @@ class GP(model):
         KiKx = np.dot(self.Ki,Kx)
         if full_cov:
             Kxx = self.kern.K(_Xnew, slices1=slices,slices2=slices)
-            var = Kxx - np.dot(KiKx.T,Kx) #NOTE is the shape of v right?
+            var = Kxx - np.dot(KiKx.T,Kx) #NOTE this won't work for plotting
         else:
             Kxx = self.kern.Kdiag(_Xnew, slices=slices)
             var = Kxx - np.sum(np.multiply(KiKx,Kx),0)
diff --git a/GPy/models/sparse_GP.py b/GPy/models/sparse_GP.py
index 6ba74e38..a90f73cb 100644
--- a/GPy/models/sparse_GP.py
+++ b/GPy/models/sparse_GP.py
@@ -196,15 +196,14 @@ class sparse_GP(GP):
 
         Kx = self.kern.K(self.Z, Xnew)
         mu = mdot(Kx.T, self.C/self.scale_factor, self.psi1V)
-
         if full_cov:
             Kxx = self.kern.K(Xnew)
-            var = Kxx - mdot(Kx.T, (self.Kmmi - self.C/self.scale_factor**2), Kx)
+            var = Kxx - mdot(Kx.T, (self.Kmmi - self.C/self.scale_factor**2), Kx) #NOTE thiswon't work for plotting
         else:
             Kxx = self.kern.Kdiag(Xnew)
             var = Kxx - np.sum(Kx*np.dot(self.Kmmi - self.C/self.scale_factor**2, Kx),0)
 
-        return mu,var
+        return mu,var[:,None]
 
     def plot(self, *args, **kwargs):
         """