)*eye(n))';
    end
    iK(:,:,i) = L'\(L\eye(n));
    beta(:,i) = L'\(L\gpmodel.targets(:,i));
  end
  
end

% initializations
k = zeros(n,E); M = zeros(E,1); V = zeros(D,E); S = zeros(E);
dMds = zeros(E,D,D); dSdm = zeros(E,E,D); r = zeros(1,D);
dSds = zeros(E,E,D,D); dVds = zeros(D,E,D,D); T = zeros(D);
tlbdi = zeros(n,D); dMdi = zeros(E,n,D); dMdt = zeros(E,n,E);
dVdt = zeros(D,E,n,E); dVdi = zeros(D,E,n,D); dSdt = zeros(E,E,n,E);
dSdi = zeros(E,E,n,D); dMdX = zeros(E,D+2,E); dSdX = zeros(E,E,D+2,E);
dV