function tvMatrixInverse(n,to,tf,j,s,eta,tau,comp)

% Inputs : n is the dimension of the test matrix (2 <= n <= 7 allowed)
%          to, tf are the time limits from to to tf
%          j, s specify the type of look-ahead convergent discretization
%               formula used
%          eta is the decay constant for the ZNN method
%          tau is the sampling gap
%          comp = 1 accuracy comparison performed of the number of correct 
%                   leading digits in every entry of computed X(t) = A^{-1}(t)

% Sample call tvMatrixInverse(7,-0.274,10,4,5,300,0.0002,1)

[polyrest, taucoeff] = Polyksrestcoeff3(j,s); 
tautau = taucoeff*tau;

pno = (tf - to)/tau;        % number of eigenanalyses for pno steps by tau
d = j + s;                  % length of polyrest for a (d+1)-IFD
dA = min(j+s,7);            % length of difference equation for Adot
pnosignif = round(pno/10,1,'significant'); % for 10 extra comparisons
kktry = [20:pnosignif:pno,pno+2];  % indices used for comparisons

Z = []; AA = [];            % initializing settings

for k = 1:d                 % compute inverses for j+s start-up matrices A(t) 
  tk = to + (k-1)*tau; 
  Atk = formA(n,tk);          %         <<== change A formation (twice) ! ! ! 
  inA = inv(Atk);             % inverses store in column form for lin combinations
  Z = [inA(:),Z];           % start-up inverse flow matrices; latest one first
  if k+dA > d
    AA = [Atk(:),AA];         % storing relevant A matrices from the front back
  end
end

for kk = d+1:pno+2                       % for each discrete time step
  Adot = reshape(dotA(n,AA,dA,tau),n,n); % Atk and Adot are n by n matrices
  Xk = reshape(Z(:,1),n,n);
  Y = -tautau*Xk*((Adot+eta*Atk)*Xk - eta*eye(n)); % this is an n by n matrix 
  W = -Z*polyrest;                       % this is in column matrix form  
  Z = [Y(:)+W,Z(:,1:end-1)];             % find new inverse 
  tk = tk + tau;                         % update time
  Atk = formA(n,tk);                     %    and A.. matrix <<== adjust ! !
  AA = [Atk(:),AA(:,1:end-1)];           %         and AA
  
  if comp == 1 && ismember(kk,kktry)     % then compare
     Xk = reshape(Z(:,1),n,n);
     format short e
     nrelerrkk = [kk,norm(Atk*Xk - eye(n))/norm(Atk)], format
  end
end
% tf = tk - tau  % check the end time tf

end






