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

% Inputs : m,n are the dimension of the test matrix A(t)_{m,n}, m >= n 
%              (1 <= m,n <= 8 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(tk)
%                   when solving  linear equation A(t)x(t) = b(t)
%                   predictively.
% External m-files used: formAmn.m, formb.m; dotAmn.m, dotb.m;
%                        Polyksrestcoeff3.m

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

if m < n 
  error(['You may be looking for a right pseudoinverse if m < n; here we work'...
          'with m > n and left pseudoinverses only'])
end

[polyrest, taucoeff] = Polyksrestcoeff3(j,s); % set up look-ahead finite 
tautau = taucoeff*tau;      %   difference formula j_s coefficients

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 13 accuracy comparisons
kktry = [d+1,20:pnosignif:pno,pno+2]; % iteration indices used in comparisons
AA = []; bb = []; X = [];   % initializing settings for data sets 

for k = 1:d                 % compute pseudoinverses for j+s start-up matrices A(tk) 
  tk = to + (k-1)*tau; 
  Atk = formAmn(m,n,tk);    %         <<== change A formation (twice) ! ! ! 
  btk = formb(m,tk);        % form b(tk)
  xk = inv(Atk'*Atk)*(Atk'*btk); % xk solves normal equation of A(tk)x(tk)=b(tk)
  X = [xk,X];               % store solutions
  if k+dA > d
    AA = [Atk(:),AA];       % store A matrices from the front on back
    bb = [btk,bb];          % store right-hand sides b(tk)
  end
end

for kk = d+1:pno+2             % for each discrete time step
  Adotmn = reshape(dotAmn(m,n,AA,dA,tau),m,n); % Atk and Adotmn are m by n matrices  
  bdot = dotb(bb,dA,tau);                    % bdot and b are m vectors
  
% xk = tautau*((Atk'*Atk)\(Atk'*(-(Adotmn + eta*Atk)*xk + bdot + eta*btk))) - X*polyrest;
%           taking the inverse sems always quicker here !
 
  xk = tautau*inv(Atk'*Atk)*Atk'*(-(Adotmn + eta*Atk)*xk + bdot + eta*btk) - X*polyrest;

  X = [xk,X(:,1:end-1)];       % apply ZNN formula for j_s formula
  tk = tk + tau;               % update time
  Atk = formAmn(m,n,tk);       %  and A.. matrix          <<== adjust ! !
  AA = [Atk(:),AA(:,1:end-1)]; %    and store Atk matrices for Atk derivative
  btk = formb(m,tk);           % update b(tk) 
  bb = [btk,bb(:,1:end-1)];    % store b(tk)
   
  if comp == 1 && ismember(kk,kktry)     % then compare if desired
    format short e
    npslinequat = [kk,norm(Atk*xk - btk)/norm(Atk)], % check solution
    format
  end
end
%  tf = tk - tau               % check that the end time tf has been reached

end
