function [ar,br,cr,dr,er] = elimAX(a,b,c,d,e,iax)
%ELIMAX  Eliminates algebraic states in DAE
%
%        E dx1/dt  = A11 x1 + A12 x2 + B1 u
%                0 = A21 x1 + A22 x2 + B2 u
%                y =  C1 x1 +  C2 x2 +  D u
%
%  IAX contains the indices of the algebraic states x2=x(IAX).
%
%  Note: Returns Dr=NaN(size(D)) when A22 is numerically singular.

%   Copyright 1986-2020 The MathWorks, Inc. 
M = [a b;c d];
nx = size(a,1);
ikeep = 1:nx;  ikeep(iax) = [];
jkeep = 1:nx;  jkeep(iax) = [];
er = e(ikeep,jkeep);
% Schur complement
[nr,nc] = size(M);
ikeep = [ikeep nx+1:nr];
jkeep = [jkeep nx+1:nc];
Mr = M(ikeep,jkeep) - M(ikeep,iax) * ...
   matlab.internal.math.nowarn.mldivide(M(iax,iax),M(iax,jkeep));

if allfinite(Mr)
   nx = nx-numel(iax);
   [nr,nc] = size(Mr);
   ar = Mr(1:nx,1:nx);
   br = Mr(1:nx,nx+1:nc);
   cr = Mr(nx+1:nr,1:nx);
   dr = Mr(nx+1:nr,nx+1:nc);
   if isequal(er,eye(nx))
      er = [];
   end
else
   [ny,nu] = size(d);
   dr = NaN(ny,nu);
   ar = [];  er = [];  br = zeros(0,nu);  cr = zeros(ny,0);
end