diff --git a/GPy/examples/state_space.py b/GPy/examples/state_space.py
index ce506b9c..6d96a927 100644
--- a/GPy/examples/state_space.py
+++ b/GPy/examples/state_space.py
@@ -7,36 +7,36 @@ import GPy.models.state_space_new as SS_new
 #X = np.linspace(0, 10, 2000)[:, None]
 #Y = np.sin(X) + np.random.randn(*X.shape)*0.1
 
-## Need to run these lines when X and Y are imported ->
-#X.shape = (X.shape[0],1)
-#Y.shape = (Y.shape[0],1)
-## Need to run these lines when X and Y are imported <-
+# Need to run these lines when X and Y are imported ->
+X.shape = (X.shape[0],1)
+Y.shape = (Y.shape[0],1)
+# Need to run these lines when X and Y are imported <-
 
-# Generation of minimal example data ->
-X = np.random.rand(3)
-sort_index = np.argsort(X)
-X = X[sort_index]; X.shape = (X.shape[0],1)
-Y = np.sin(10*X) + np.random.randn(*X.shape)*0.1
-# Generation of minimal example data <-
+## Generation of minimal example data ->
+#X = np.random.rand(3)
+#sort_index = np.argsort(X)
+#X = X[sort_index]; X.shape = (X.shape[0],1)
+#Y = np.sin(10*X) + np.random.randn(*X.shape)*0.1
+## Generation of minimal example data <-
 
 #plt.figure()
 #plt.plot( X, Y)
 #plt.show()
 
-kernel = GPy.kern.Matern32(X.shape[1])
-m = GPy.models.StateSpace(X,Y, kernel)
-
-print m
+#kernel = GPy.kern.Matern32(X.shape[1])
+#m = GPy.models.StateSpace(X,Y, kernel)
 #
-m.optimize()
-#
-print m
+#print m
+##
+#m.optimize(optimizer='bfgs',messages=True)
+##
+#print m
 
 kernel1 = GPy.kern.Matern32(X.shape[1])
 m1  = GPy.models.GPRegression(X,Y, kernel1)
 
 print m1
-m1.optimize()
+m1.optimize(optimizer='bfgs',messages=True)
 
 print m1
 
@@ -45,7 +45,7 @@ m2  = SS_new.StateSpace(X,Y, kernel2)
 
 print m2
 
-m2.optimize()
+m2.optimize(optimizer='bfgs',messages=True)
 
 print m2
 
diff --git a/GPy/models/state_space_main.py b/GPy/models/state_space_main.py
index d6260bc3..b903ef34 100644
--- a/GPy/models/state_space_main.py
+++ b/GPy/models/state_space_main.py
@@ -657,6 +657,7 @@ class DescreteStateSpace(object):
                 dP_pred[:,:,j] += dP_pred[:,:,j].T            
                 dP_pred[:,:,j] += np.dot( A ,np.dot(dP, A.T)) + dQ
                 
+                dP_pred[:,:,j] = 0.5*(dP_pred[:,:,j] + dP_pred[:,:,j].T) #symmetrize
         else:
             dm_pred = None
             dP_pred = None
@@ -1439,7 +1440,7 @@ class ContDescrStateSpace(DescreteStateSpace):
         unique_round_decimals = 8
         dt = np.empty((X.shape[0],))
         dt[1:] = np.diff(X[:,0],axis=0)
-        dt[0]  = dt[1]
+        dt[0]  = 0#dt[1]
         unique_indices = np.unique(np.round(dt, decimals=unique_round_decimals))
     
         if len(unique_indices) > 20:        
@@ -1531,13 +1532,19 @@ class ContDescrStateSpace(DescreteStateSpace):
  
             # The dynamical model
             A  = linalg.expm(F*dt) 
- 
+            if np.any( np.isnan(A)):
+                A  = linalg.expm3(F*dt) 
+            
             # The covariance matrix Q by matrix fraction decomposition ->
             Phi = np.zeros((2*n,2*n))
             Phi[:n,:n] = F
             Phi[:n,n:] = L.dot(Qc).dot(L.T)
             Phi[n:,n:] = -F.T
-            AB = linalg.expm(Phi*dt).dot(np.vstack((np.zeros((n,n)),np.eye(n))))
+            AB = linalg.expm(Phi*dt)
+            if np.any(np.isnan(AB)):
+                AB = linalg.expm3(Phi*dt)
+            AB = np.dot(AB, np.vstack((np.zeros((n,n)),np.eye(n))))
+            
             Q_noise_1 = linalg.solve(AB[n:,:].T,AB[:n,:].T)
             # The covariance matrix Q by matrix fraction decomposition <-