inp2(:,d)*inp2(:,d)';
    dKdX = -K(:,:,i).*dsqdX./2;                               <span class="comment">% dK/dX(1:D)</span>
    dKdXbeta = dKdX*beta(:,i);
    bdX(:,i,d) = -K2\dKdXbeta;                        <span class="comment">% dbeta/dX</span>
    dslb(d) = -liK'*dKdXbeta + beta(:,i)'*dldX(:,d);
    dlb = dldX(:,d).*beta(:,i) + l.*bdX(:,i,d);
    dtdX = inp*(-bsxfun(@times,iR(:,d),2.*exp(2*X(d,i))*iR(d,:)));
    dlbt = lb'*dtdX + dlb'*t;
    dVdX(:,i,d,i) = (dlbt'*c + cdX(d)*(lb'*t)');
  <span class=