def sde(self):
        """
        Return the state space representation of the covariance.
        
        Note! For Sparse GP inference too small or two high values of lengthscale
        lead to instabilities. This is because Qc are too high or too low
        and P_inf are not full rank. This effect depends on approximatio order.
        For N = 10. lengthscale must be in (0.8,8). For other N tests must be conducted.
        N=6: (0.06,31)
        Variance should be within reasonable bounds as well, but its dependence is linear.
        
        The above facts do not take into accout regularization.
        """
        #import pdb; pdb.set_trace()
        if self.approx_order is not None:
            N = self.approx_order
        else:
            N = 10# approximation order ( number of terms in exponent series expansion)
            
        roots_rounding_decimals = 6

        fn = np.math.factorial(N)

        p_lengthscale = float( self.lengthscale )
        p_variance = float(self.variance)
        kappa = 1.0/2.0/p_lengthscale**2

        Qc = np.array( ((p_variance*np.sqrt(np.pi/kappa)*fn*(4*kappa)**N,),) )
        
        eps = 1e-12
        if (float(Qc) > 1.0/eps) or (float(Qc) < eps):
            warnings.warn("""sde_RBF kernel: the noise variance Qc is either very large or very small. 
                                It influece conditioning of P_inf: {0:e}""".format(float(Qc)) )

        pp1 = np.zeros((2*N+1,)) # array of polynomial coefficients from higher power to lower

        for n in range(0, N+1): # (2N+1) - number of polynomial coefficients
            pp1[2*(N-n)] = fn*(4.0*kappa)**(N-n)/np.math.factorial(n)*(-1)**n
            
        pp = sp.poly1d(pp1)
        roots = sp.roots(pp)

        neg_real_part_roots = roots[np.round(np.real(roots) ,roots_rounding_decimals) < 0]
        aa = sp.poly1d(neg_real_part_roots, r=True).coeffs

        F = np.diag(np.ones((N-1,)),1)
        F[-1,:] = -aa[-1:0:-1]

        L= np.zeros((N,1))
        L[N-1,0] = 1

        H = np.zeros((1,N))
        H[0,0] = 1

        # Infinite covariance:
        Pinf = lyap(F, -np.dot(L,np.dot( Qc[0,0],L.T)))
        Pinf = 0.5*(Pinf + Pinf.T)
        # Allocating space for derivatives
        dF    = np.empty([F.shape[0],F.shape[1],2])
        dQc   = np.empty([Qc.shape[0],Qc.shape[1],2])
        dPinf = np.empty([Pinf.shape[0],Pinf.shape[1],2])

        # Derivatives:
        dFvariance = np.zeros(F.shape)
        dFlengthscale = np.zeros(F.shape)
        dFlengthscale[-1,:] = -aa[-1:0:-1]/p_lengthscale * np.arange(-N,0,1)

        dQcvariance = Qc/p_variance
        dQclengthscale = np.array(( (p_variance*np.sqrt(2*np.pi)*fn*2**N*p_lengthscale**(-2*N)*(1-2*N),),))
        
        dPinf_variance = Pinf/p_variance

        lp = Pinf.shape[0]
        coeff = np.arange(1,lp+1).reshape(lp,1) + np.arange(1,lp+1).reshape(1,lp) - 2
        coeff[np.mod(coeff,2) != 0] = 0
        dPinf_lengthscale = -1/p_lengthscale*Pinf*coeff

        dF[:,:,0]    = dFvariance
        dF[:,:,1]    = dFlengthscale
        dQc[:,:,0]   = dQcvariance
        dQc[:,:,1]   = dQclengthscale
        dPinf[:,:,0] = dPinf_variance
        dPinf[:,:,1] = dPinf_lengthscale

        P0 = Pinf.copy()
        dP0 = dPinf.copy()

        if self.balance:
            # Benefits of this are not very sound. Helps only in one case:
            # SVD Kalman + RBF kernel
            import GPy.models.state_space_main as ssm
            (F, L, Qc, H, Pinf, P0, dF, dQc, dPinf,dP0) = ssm.balance_ss_model(F, L, Qc, H, Pinf, P0, dF, dQc, dPinf, dP0 )

        return (F, L, Qc, H, Pinf, P0, dF, dQc, dPinf, dP0)