diff --git a/GPy/kern/parts/ODE_UY.py b/GPy/kern/parts/ODE_UY.py
index bb736cc5..3ddf174b 100644
--- a/GPy/kern/parts/ODE_UY.py
+++ b/GPy/kern/parts/ODE_UY.py
@@ -186,20 +186,29 @@ class ODE_UY(Kernpart):
 
     def dK_dtheta(self, dL_dK, X, X2, target):
         """derivative of the covariance matrix with respect to the parameters."""
-        if X2 is None: X2 = X
-        dist = np.abs(X - X2.T)
-
+       
         X,slices = X[:,:-1],index_to_slices(X[:,-1])
         if X2 is None:
             X2,slices2 = X,slices
         else:
             X2,slices2 = X2[:,:-1],index_to_slices(X2[:,-1])
-
-
+        #rdist = X[:,0][:,None] - X2[:,0][:,None].T
+        rdist = X - X2.T
         ly=1/self.lengthscaleY
         lu=np.sqrt(3)/self.lengthscaleU
-        #ly=self.lengthscaleY
-        #lu=self.lengthscaleU
+
+        rd=rdist.shape[0]
+        dktheta1 = np.zeros([rd,rd])
+        dktheta2 = np.zeros([rd,rd])
+        dkdvar = np.zeros([rd,rd])
+
+        # dk dtheta for UU
+        UUdtheta1 = lambda dist: np.exp(-lu* dist)*dist + (-dist)*np.exp(-lu* dist)*(1+lu*dist)
+        UUdtheta2 = lambda dist: 0
+        UUdvar = lambda dist: (1 + lu *dist) * np.exp(-lu* dist)
+
+
+        # dk dtheta for YY
 
         dk1theta1 = lambda dist: np.exp(-ly*dist)*2*(-lu)/(lu+ly)**3
         #c=np.sqrt(3)
@@ -216,7 +225,7 @@ class ODE_UY(Kernpart):
       
         dk3theta1 = lambda dist: np.exp(-dist*lu)*(lu+ly)**(-2)*((2*lu+ly+dist*lu**2+lu*ly*dist)*(-dist-2/(lu+ly))+2+2*lu*dist+ly*dist)
 
-        dktheta1 = lambda dist: self.varianceU*self.varianceY*(dk1theta1+dk2theta1+dk3theta1)
+        #dktheta1 = lambda dist: self.varianceU*self.varianceY*(dk1theta1+dk2theta1+dk3theta1)
 
 
 
@@ -230,14 +239,20 @@ class ODE_UY(Kernpart):
 
         dk3theta2 = lambda dist: np.exp(-dist*lu) * (-3*lu-ly-dist*lu**2-lu*ly*dist)/(lu+ly)**3
 
-        dktheta2 = lambda dist: self.varianceU*self.varianceY*(dk1theta2 + dk2theta2 +dk3theta2)
+        #dktheta2 = lambda dist: self.varianceU*self.varianceY*(dk1theta2 + dk2theta2 +dk3theta2)
 
 
 
         k1 = lambda dist: np.exp(-ly*dist)*(2*lu+ly)/(lu+ly)**2
         k2 = lambda dist: (np.exp(-lu*dist)*(ly-2*lu+lu*ly*dist-lu**2*dist) + np.exp(-ly*dist)*(2*lu-ly) ) / (ly-lu)**2 
         k3 = lambda dist: np.exp(-lu*dist) * ( (1+lu*dist)/(lu+ly) + (lu)/(lu+ly)**2 )
-        dkdvar = k1+k2+k3
+        #dkdvar = k1+k2+k3
+        
+
+        # dk dtheta for UY
+
+
+
 
 
         for i, s1 in enumerate(slices):
@@ -246,16 +261,35 @@ class ODE_UY(Kernpart):
                     for ss2 in s2:
                         if i==0 and j==0:
                             #target[ss1,ss2] = kuu(np.abs(rdist[ss1,ss2]))
+                            #dktheta1[ss1,ss2] =
+                            #dktheta2[ss1,ss2] =
+                            #dkdvar[ss1,ss2] =
+                            dktheta1[ss1,ss2] = self.varianceU*self.varianceY*UUdtheta1(rdist[ss1,ss2])
+                            dktheta2[ss1,ss2] = 0
+                            dkdvar[ss1,ss2] = self.varianceY*UUdvar(rdist[ss1,ss2])
                         elif i==0 and j==1:
                             #target[ss1,ss2] = np.where(  rdist[ss1,ss2]>0 , kuyp(np.abs(rdist[ss1,ss2])), kuyn(np.abs(rdist[s1[0],s2[0]]) )   )
+                            #dktheta1[ss1,ss2] =
+                            #dktheta2[ss1,ss2] =
+                            #dkdvar[ss1,ss2] =                            
+                            dktheta1[ss1,ss2] = self.varianceU*self.varianceY*(dk1theta1(np.abs(rdist[ss1,ss2]))+dk2theta1(np.abs(rdist[ss1,ss2]))+dk3theta1(np.abs(rdist[ss1,ss2])))
+                            dktheta2[ss1,ss2] = self.varianceU*self.varianceY*(dk1theta2(np.abs(rdist[ss1,ss2])) + dk2theta2(np.abs(rdist[ss1,ss2])) +dk3theta2(np.abs(rdist[ss1,ss2])))
+                            dkdvar[ss1,ss2] = k1(np.abs(rdist[ss1,ss2]))+k2(np.abs(rdist[ss1,ss2]))+k3(np.abs(rdist[ss1,ss2]))
                         elif i==1 and j==1:
                             #target[ss1,ss2] = kyy(np.abs(rdist[ss1,ss2]))
+                            dktheta1[ss1,ss2] = self.varianceU*self.varianceY*(dk1theta1(np.abs(rdist[ss1,ss2]))+dk2theta1(np.abs(rdist[ss1,ss2]))+dk3theta1(np.abs(rdist[ss1,ss2])))
+                            dktheta2[ss1,ss2] = self.varianceU*self.varianceY*(dk1theta2(np.abs(rdist[ss1,ss2])) + dk2theta2(np.abs(rdist[ss1,ss2])) +dk3theta2(np.abs(rdist[ss1,ss2])))
+                            dkdvar[ss1,ss2] = k1(np.abs(rdist[ss1,ss2]))+k2(np.abs(rdist[ss1,ss2]))+k3(np.abs(rdist[ss1,ss2]))
                         else:
                             #target[ss1,ss2] = np.where(  rdist[ss1,ss2]>0 , kyup(np.abs(rdist[ss1,ss2])), kyun(np.abs(rdist[s1[0],s2[0]]) )   )
+                            #dktheta1[ss1,ss2] =
+                            #dktheta2[ss1,ss2] =
+                            #dkdvar[ss1,ss2] =                            
+                            dktheta1[ss1,ss2] = self.varianceU*self.varianceY*(dk1theta1(np.abs(rdist[ss1,ss2]))+dk2theta1(np.abs(rdist[ss1,ss2]))+dk3theta1(np.abs(rdist[ss1,ss2])))
+                            dktheta2[ss1,ss2] = self.varianceU*self.varianceY*(dk1theta2(np.abs(rdist[ss1,ss2])) + dk2theta2(np.abs(rdist[ss1,ss2])) +dk3theta2(np.abs(rdist[ss1,ss2])))
+                            dkdvar[ss1,ss2] = k1(np.abs(rdist[ss1,ss2]))+k2(np.abs(rdist[ss1,ss2]))+k3(np.abs(rdist[ss1,ss2]))
 
-
-
-
+        #stop                   
 
 
         target[0] += np.sum(self.varianceY*dkdvar * dL_dK)