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

% Inputs : n is the dimension of the test matrix $A(t) (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)
%                   and in every linear equation's solution
% External m-files used: formA.m, formb.m; dotA.m, dotb.m;
%                        Polyksrestcoeff3.m

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

[polyrest, taucoeff] = Polyksrestcoeff3(j,s); % get look ahead finite 
tautau = taucoeff*tau;      %             difference  formyula coeficients

pno = (tf - to)/tau;        % number of steps with tau from to to tf
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 data sets
bb = []; Y = [];  

for k = 1:d                 % compute inverses for j+s start-up matrices A(tk) 
  tk = to + (k-1)*tau; 
  Atk = formA(n,tk);        %         <<== change A formation (twice) ! ! ! 
  btk = formb(n,tk);        %         <<== change b formation (twice) ! ! ! 
  yk = Atk\btk;
  Y = [yk,Y];               % solution vectors to linear equation A(tk) x(tk) = b(tk)
  inA = inv(Atk);           % matr inverses (as columns for forming 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
    bb = [btk,bb];          % right-hand side vectors from newest on 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 
  bdot = dotb(bb,dA,tau);              % derivatives of both at time tk
  Xk = reshape(Z(:,1),n,n);              % matrix of inverse of A(tk)

  yk = tautau*(Xk*(-Adot*yk+bdot+eta*btk)-eta*yk)-Y*polyrest; % l eq solution 
  Y = [yk,Y(:,1:end-1)];                 % future lin equat solution at tk + tau

  V = -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 = [V(:)+W,Z(:,1:end-1)];             % new inverse of A(tk+1)
 
  tk = tk + tau;                         % update time
  Atk = formA(n,tk);                     %    and A.. matrix <<== adjust ! !
  AA = [Atk(:),AA(:,1:end-1)];           %         and AA
  btk = formb(n,tk);                     % same for b of lin equation
  bb = [btk,bb(:,1:end-1)];
  
  if comp == 1 && ismember(kk,kktry)     % then compare if desired
    Xk = reshape(Z(:,1),n,n);
    format short e
    nrelerrkkXkAxb = [kk,norm(Atk*Xk - eye(n))/norm(Atk),...
         norm(Atk*yk - btk)/norm(Atk)],
    format
  end
end
% tf = tk - tau               % check that the end time tf has been reached

end






