diff --git a/GPy/examples/regression.py b/GPy/examples/regression.py
index 7d092c26..90fca886 100644
--- a/GPy/examples/regression.py
+++ b/GPy/examples/regression.py
@@ -123,7 +123,7 @@ def coregionalisation_toy():
     Y = np.vstack((Y1,Y2))
 
     k1 = GPy.kern.rbf(1)
-    k2 = GPy.kern.coregionalise(2,1)
+    k2 = GPy.kern.coregionalise(2,2)
     k = k1.prod_orthogonal(k2)
     m = GPy.models.GP_regression(X,Y,kernel=k)
     m.constrain_fixed('rbf_var',1.)
@@ -147,8 +147,8 @@ def coregionalisation_sparse():
     """
     A simple demonstration of coregionalisation on two sinusoidal functions
     """
-    X1 = np.random.rand(500,1)*8
-    X2 = np.random.rand(300,1)*5
+    X1 = np.random.rand(50,1)*8
+    X2 = np.random.rand(30,1)*5
     index = np.vstack((np.zeros_like(X1),np.ones_like(X2)))
     X = np.hstack((np.vstack((X1,X2)),index))
     Y1 = np.sin(X1) + np.random.randn(*X1.shape)*0.05
@@ -158,7 +158,7 @@ def coregionalisation_sparse():
     Z = np.hstack((np.random.rand(25,1)*8,np.random.randint(0,2,25)[:,None]))
 
     k1 = GPy.kern.rbf(1)
-    k2 = GPy.kern.coregionalise(2,2)
+    k2 = GPy.kern.coregionalise(2,1)
     k = k1.prod_orthogonal(k2) + GPy.kern.white(2,0.001)
 
     m = GPy.models.sparse_GP_regression(X,Y,kernel=k,Z=Z)
@@ -180,7 +180,6 @@ def coregionalisation_sparse():
     return m
 
 
-
 def multiple_optima(gene_number=937,resolution=80, model_restarts=10, seed=10000):
     """Show an example of a multimodal error surface for Gaussian process regression. Gene 939 has bimodal behaviour where the noisey mode is higher."""
 
diff --git a/GPy/kern/coregionalise.py b/GPy/kern/coregionalise.py
index c24cb568..f6b9426f 100644
--- a/GPy/kern/coregionalise.py
+++ b/GPy/kern/coregionalise.py
@@ -4,6 +4,7 @@
 from kernpart import kernpart
 import numpy as np
 from GPy.util.linalg import mdot, pdinv
+import pdb
 
 class coregionalise(kernpart):
     """
@@ -37,7 +38,6 @@ class coregionalise(kernpart):
         self.B = np.dot(self.W,self.W.T) + np.diag(self.kappa)
 
     def _get_param_names(self):
-
         return sum([['W%i_%i'%(i,j) for j in range(self.R)] for i in range(self.Nout)],[]) + ['kappa_%i'%i for i in range(self.Nout)]
 
     def K(self,index,index2,target):
@@ -46,7 +46,8 @@ class coregionalise(kernpart):
             index2 = index
         else:
             index2 = np.asarray(index2,dtype=np.int)
-        ii,jj = np.meshgrid(index2,index)
+        ii,jj = np.meshgrid(index,index2)
+        ii,jj = ii.T, jj.T
         target += self.B[ii,jj]
 
     def Kdiag(self,index,target):
@@ -58,19 +59,22 @@ class coregionalise(kernpart):
             index2 = index
         else:
             index2 = np.asarray(index2,dtype=np.int)
-        ii,jj = np.meshgrid(index2,index)
-        PK = np.zeros((self.R,self.R))
+        ii,jj = np.meshgrid(index,index2)
+        ii,jj = ii.T, jj.T
+
         partial_small = np.zeros_like(self.B)
         for i in range(self.Nout):
-            for j in range(self.Nout):
-                partial_small[i,j] = np.sum(partial[(ii==i)*(jj==j)])
-        dkappa = np.diag(partial_small)
+            for j in range(i,self.Nout):
+                tmp = np.sum(partial[(ii==i)*(jj==j)])
+                partial_small[i,j] = tmp
+                partial_small[j,i] = tmp
 
+        dkappa = np.diag(partial_small)
         dW = 2.*(self.W[:,None,:]*partial_small[:,:,None]).sum(0)
 
         target += np.hstack([dW.flatten(),dkappa])
 
-    def dKdiag_dtheta_foo(self,partial,index,target):
+    def dKdiag_dtheta(self,partial,index,target):
         index = np.asarray(index,dtype=np.int).flatten()
         partial_small = np.zeros(self.Nout)
         for i in range(self.Nout):
@@ -82,7 +86,5 @@ class coregionalise(kernpart):
     def dK_dX(self,partial,X,X2,target):
         pass
 
-    def dKdiag_dtheta(self,partial,index,target):
-        self.dK_dtheta(np.diag(partial),index,index,target)