) = beta(:,i)'*L/(K(:,:,j)+exp(2*X(D+2,j))*eye(n))*t <span class="keyword">...</span>
        - dMdt(j,:,j)*M(i);
      dSdt(j,i,:,:) = dSdt(i,j,:,:);
      dSdX(i,j,:,j) = reshape(dSdX(i,j,:,j),1,D+2) - M(i)*dMdX(j,:,j);
      dSdX(i,j,:,i) = reshape(dSdX(i,j,:,i),1,D+2) - M(j)*dMdX(i,:,i);
    <span class="keyword">end</span>

    dSdi(i,j,:,:) = Z*t - reshape(M(i)*dMdi(j,:,:) + dMdi(i,:,:)*M(j),n,D);
    dSdi(j,i,:,:) = dSdi(i,j,:,:);
    dSdX(j,i,:,:) = dSdX(i,j,:,:);
  <span class="keyword">end</span> 