p(:,d)');
      B = K(:,:,i).*Q;
      Z(:,d) = Z(:,d)+exp(-2*X(d,i))*(B*beta(:,i).*bLiKj+beta(:,i).*(B*bLiKj));

      <span class="keyword">if</span> i~=j; B = K(:,:,j).*Q; <span class="keyword">end</span>

      Z(:,d) = Z(:,d)+exp(-2*X(d,j))*(bLiKi.*(B*beta(:,j))+B*bLiKi.*beta(:,j));
      B = bsxfun(@plus,zi(:,d),zj(:,d)').*A;
      r(d) = sum(sum(B))*t;
      T(d,1:d) = sum(zi(:,1:d)'*B,2) + sum(B*zj(:,1:d))';
      T(1:d,d) = T(d,1:d)';

      <span class="keyword">if</span> i==j
        RTi =  bsxfu