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

% Inputs : m,n are the dimension of the test matrix A(t)_{m,n} (2 <= 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(t) = A^{-1}(t)
%                   and in every linear equation's solution
% External m-files used: formAmn.m, formb.m; dotA.m, dotb.m;
%                        Polyksrestcoeff3.m

% Sample call tvLeftPseudInv(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);    % get look-ahead finite 
tautau = taucoeff*tau;      %             difference formula 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 10 extra comparisons
kktry = [d+1,20:pnosignif:pno,pno+2];  % indices used for comparisons
AA = []; G = [];            % 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) ! ! ! 
  gk = inv(Atk'*Atk)*Atk';  % matrix leftt pseudoinverse 
  G = [gk(:),G];            % store pseudo inverses from front to back
  if k+dA > d
    AA = [Atk(:),AA];       % store A matrices from the front on back
  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  
  Y = tautau*inv(Atk'*Atk)*(-((Adotmn'+eta*Atk')*Atk+Atk'*Adotmn)*gk+Adotmn'+eta*Atk');
  
%  Y = tautau*((Atk'*Atk)\(-((Adotmn'+eta*Atk')*Atk+Atk'*Adotmn)*gk+Adotmn'+eta*Atk'));
%        using inv instead of \ is quicker ! 

  W = -G*polyrest;
  G = [Y(:)+W,G(:,1:end-1)];   % construct and store sequence of pseudoinverses

  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
  gk = reshape(G(:,1),n,m);    % G is the n by m  pseudoinverse Atk^+ of Atk
  
  if comp == 1 && ismember(kk,kktry)     % then compare if desired
    format short e
    nleftpsinv = [kk,norm(gk - inv(Atk'*Atk)*Atk')/norm(Atk)],
    format
  end
end
%  tf = tk - tau               % check that the end time tf has been reached

end
