function  tvMatrSquareRootwEulerStartv(n,to,tf,j,s,eta,tau,ploterror)

% Sample call :  tvMatrSquareRootwEulerStartv(5,0,360,4,5,2,0.02,1)
%                    for relative errors of magnitude 10^-12 after 1 minute
%                    and of magnitude 10^-13 after 6 minutes

% inputs : n the dimension of the input square matrix time-varying data flow
%          [to, tf] : the time interval under consideration
%          j, s : defining integers for desired look-ahead convergent fin. diff. scheme
%          eta : the decay constant chosen for the error function DE
%          tau : the constant sampling step in digital fractions of a
%                second of time; (50 Hz = 0.02 sec for standard sensor input)
%          ploterror = 1 : this will compute the difference between A(tk) and the 
%                          computed square root X(tk) squared for each time instance
%                   ~= 1 : no errror computations
% required subcodes: A function code or sensor iput for A(tk) and tk in [to, tf].
%                    We use  formAnnv.m  here for dimensions 1 <= n <= 8.
%                    And dotAv.m for the backward finite difference formulation 
%                    of \dot Av, as well as  Polyksrestcoeff3.m  for the look-ahead 
%                    finite difference formula coefficients
% 

if nargin == 7, ploterror = 0; end            % default
warning off
eta0 = 0.0375/tau; % h ~ 0.0375 for starting values with Euclid \dot x = (x(2)-x(1))/tau

[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 backward difference equation for Adv
In = eye(n);
AA = []; X = []; nev = [];  % initializing settings for data sets 
                            % X and AA store newest entry in front, going back in time;
                            % nev stores the error value conformally with time flow

      % Starting values generation from random start and simple Euler ZNN steps
      
tk = to;                    % first two start-up values for A and X with A \approx X X
Av = formAnnv(n,tk); 
Am = reshape(Av,n,n); 
%Xm = sqrtm(Am);            % try this near-perfect first solution Xm(to) for A(to) also and see
Xm = randn(n,n);            % total random entry matrix as first 'solution' starting value
%Xm = Am.^(0.7);            % .^0.6  or .^0.7 work 'well' as 'entry' matrix Am roots for Xm and
Xv = Xm(:);                 %  the matr square root problem, but not better than a random start.        
Xsv = kron(In,Xm)*Xv; 
if ploterror == 1, nev = [nev,[norm(Av - Xsv,2)/norm(Av,2); tk]]; end
X = [Xv,X];
AA = [Av,AA];

tk = tk + tau;
Av = formAnnv(n,tk);
Adv = (Av - AA(:,1))/tau;     
Xv = X(:,1) + tau*((kron(Xm.',In) + kron(In,Xm))\(Adv + eta0*AA(:,1) - eta0*Xsv));
X = [Xv,X];
AA = [Av,AA];
Xm = reshape(Xv,n,n);        Xsv = kron(In,Xm)*Xv;
if ploterror == 1, nev = [nev,[norm(Av - Xsv,2)/norm(Av,2);tk]]; end

for k = 3:d                  % compute remaining j+s-2 start-up values for X from input
  tk = tk + tau;
  Av = formAnnv(n,tk);
  Adv = (Av - AA(:,1))/tau; 
  Xv = X(:,1) + tau*((kron(Xm.',In) + kron(In,Xm))\(Adv + eta0*AA(:,1) - eta0*Xsv));
  X = [Xv,X]; 
  AA = [Av,AA];
  Xm = reshape(Xv,n,n);        Xsv = kron(In,Xm)*Xv; 
  if ploterror == 1, nev = [nev,[norm(Av - Xsv,2)/norm(Av,2);tk]]; end      
end  

AA = AA(:,1:dA);             % shrink AA down to proper size n^2 by dA from n^2 by d
    
for kk = d+1:pno+2 % iterate with chosen look-ahead convergent finite difference formula
  tk = tk + tau;
  Av = formAnnv(n,tk); 
  Adv = dotAv(AA,dA,tau);
%  Xv = tautau*((kron(Xm.',In) + kron(In,Xm))\(Adv + eta*AA(:,1) - eta*Xsv)) - X*polyrest;
                    %  backslash \ version above is slower than linsolve(.. , ..)  ! !                        
  Xv = tautau*linsolve(kron(Xm.',In) + kron(In,Xm),Adv + eta*AA(:,1) - eta*Xsv) - X*polyrest;
  X = [Xv,X(:,1:end-1)]; 
  AA = [Av,AA(:,1:end-1)];
  Xm = reshape(Xv,n,n);        Xsv = kron(In,Xm)*Xv; 
  if ploterror == 1, nev = [nev,[norm(Av - Xsv,2)/norm(Av,2);tk]]; end   
end

if ploterror == 1             % plot error if desired
  figure(1)
  semilogy(nev(2,:),nev(1,:),'-b'), grid on,  
  title(['Rel. errors of time-varying matrix square roots; n = ',num2str(n),...
      ', finite difference formula j\_s = ', num2str(j),...
      '\_',num2str(s),', \eta = ',num2str(eta),', \tau = ',num2str(tau)],'FontSize',20)
  xlabel(['Seconds from to = ',num2str(to),' to tf = ',num2str(tf)],'FontSize',18), hold off
  hold off
end

