function Zmatrixeig3_3b(n,zeich,jump1,jump2,eta,tau) 
% For time-varying matrix eigenalanlyses via 1-step ahead discretization
% formula

% via 9-IFD difference equations for all eigenvalues, one at a time

% sample call :  Zmatrixeig3_3b(5,1, 4, 13, 12,0.001)

% n = size of square A(t) matrix, zeich = 1 means plot, else no plots
% jump1 is the time at first data jump, jump2 at second data jump back to original
% function . jump1 <= jump2 is necessary here. 
% (if jump2 < jump1, then the initial function is used throughout: no data jump.
% Both may lie outside the time interval [0, tf].

% call >> Zmatrixeig3_3b(5,1,2,7.5,40)   
% This function needs to define two different input A(t) functions as done in 
% inputeigsinglewjumps.m



if nargin == 1                                            % defaults
  zeich = 0; jump1 = 20; jump2 = 20; eta = 20; tau = 0.001;
end
if nargin == 2 
  jump1 = 20; jump2 = 20; eta = 20; tau = 0.001;
end
if nargin == 3
  jump2 = 20; eta = 20; tau = 0.001;
end
if nargin == 4
  eta = 20; tau = 0.001; 
end
if nargin == 5
   tau = 0.001;
end
warning('OFF'), 

tf = 20;             % adjust  tf to taste and problem!
if jump1 > 0, alt = 0; else, alt = 1; end

[f,fj,n] = inputeigsinglewjumps2(n);  % prepare input functions f = A(t) externally
if zeich == 1
  VcOrthErr = NaN(ceil(tf/tau),2);
  relerrinddat = NaN(ceil(tf/tau),n+1);
  errdat = NaN(ceil(tf/tau),2);   % prepare storage of x and lambda iterates
  evaluedat = NaN(ceil(tf/tau),n+1);
  nAdot = NaN(ceil(tf/tau),3);
else
  relerrinddat = []; 
  errdat = [];
  evaluedat = [];
end

% starting to solve the eigenproblems A(t) x(t) - lambda(t) x(t) = 0    
fprintf('     Timing for difference formulas alg (Frank U): \n')

tic % starting values for eigenvectors and eigenvalues via Francis QR :

[A1,A2,A3,A4,A5,A6,Z1,Z2,Z3,Z4,Z5,Z6,jj,relerrinddat,errdat,evaluedat] = ...
              starting(1,alt,tau,zeich,relerrinddat,errdat,evaluedat);

while jj < ceil(tf/tau)-1  % working for all time instances until tf
  jjtau = jj*tau; 
%  Adot = (147*A8-360*A7+450*A6-400*A5+225*A4-72*A3+10*A2)/(60*tau); 
  Adot = (137*A6-300*A5+300*A4-200*A3+75*A2-12*A1)/(60*tau); 
%  Adot = (25*A6-48*A5+36*A4-16*A3+3*A2)/(12*tau); 
  normAdot = norm(Adot);  
  if zeich == 1, nAdot(jj,:,:) = [jjtau,normAdot,alt]; end
 
  if normAdot > 300 % restart with four initial values for jumped function
    if jjtau > jump1 && jjtau <= jump2 % decide which function A(t) to use
      alt = 1; else, alt = 0; 
    end
    [A1,A2,A3,A4,A5,A6,Z1,Z2,Z3,Z4,Z5,Z6,jj,relerrinddat,errdat,evaluedat] = ...
                      starting(jj,alt,tau,zeich,relerrinddat,errdat,evaluedat); 
%    Adot = (147*A8-360*A7+450*A6-400*A5+225*A4-72*A3+10*A2)/(60*tau);  % use new Adot after jump
    Adot = (137*A6-300*A5+300*A4-200*A3+75*A2-12*A1)/(60*tau);
%    Adot = (25*A6-48*A5+36*A4-16*A3+3*A2)/(12*tau);
    jjtau = jj*tau;
    if zeich == 1, nAdot(jj,:,:) = [jjtau,normAdot,alt]; end 
  end   
  
  for k = 1:n   % going through each eigendata calculations for one fixed time
    z1 = Z1(:,k); z2 = Z2(:,k); z3 = Z3(:,k); z4 = Z4(:,k); z5 = Z5(:,k); 
          z6 = Z6(:,k);  % all are n+1 vectors
    z6s = z6(1:n);      % most recent kth eigenvector
    z6e = z6(n+1);      % most recent kth eigenvalue in row n+1
    
    Al = A6 - z6e*eye(n); %  sAl= size(Al)       
    P = [Al, -z6s; 2*z6s',0]; %sP = size(P)   % using Gaussian elimination with matrix P in formula (7) p.5

    q = -588*tau*[(eta*Al + Adot)*z6s; eta*(z6s'*z6s-1)]; % n+1 vector q in (7) p.5
   
    z7 = ((P\q) - 80*z6 + 182*z5 + 206*z4 - z3 - 110*z2 + 40*z1)/237;
    
    Z7(:,k) = z7;
  end           % end local eigenvalue # k loop
  
  Z1 = Z2; Z2 = Z3; Z3 = Z4; Z4 = Z5; Z5 = Z6; Z6 = Z7;  
  A1 = A2; A2 = A3; A3 = A4; A4 = A5; A5 = A6;    % cycle indices back on A1, ..., A8
  if jjtau > jump1 && jjtau <= jump2  % with discontinuous A function: 
      % here with two jumps: after t = jump1 into fj; and after t = jump2 back to f
    A6 = fj(jjtau);    % use jump function
  else
    A6 = f(jjtau);     % use standard function
  end
  
  if zeich == 1
    evaluedat(jj+1,:) = [jjtau,Z4(n+1,:)]; Vc = Z4(1:n,:);
    VcOrthErr(jj,:) = [jjtau,norm(Vc'*Vc - eye(n))];
    for k1 = 1:n
      relerrinddat(jj+1,[1,k1+1]) = [jjtau,norm(A4*Vc(:,k1) - Vc(:,k1)*Z4(n+1,k1))/abs(Z4(n+1,k1))];
    end
    errdat(jj+1,:) = [jjtau,norm(A4*Vc - Vc*diag(evaluedat(jj+1,2:end)))/norm(A4)];
  end
  jj = jj + 1;
end   % end global while loop

toc     % end of solving the  A(t) x(t) - lambda(t) x(t) eigenvalue equations

  fprintf('\n')
 
if zeich == 1
  figure(1)       % plot all n eigenvalues of A(t) for all t
  for jj = 1:n
    plot(evaluedat(:,1),evaluedat(:,jj+1),'LineWidth',0.7, 'Displayname', ['eigenvalue ',num2str(jj)])
    hold on, grid on
  end  
  a = ylim; plot([jump1,jump1],a,'k','LineWidth',1,'Displayname', ['jump1 at ',num2str(jump1)]), 
  plot([jump2,jump2],a,'k','LineWidth',1,'Displayname', ['jump2 at ',num2str(jump2)] ),
  hold off, 
  title(['Eigenvalues via 7ifd difference equations, jumps at t = ',num2str(jump1),', ',num2str(jump2)]), 
  xlabel([' 0 \leq t \leq ',num2str(tf),',  \eta = ',num2str(eta),',  \tau = ',num2str(tau)])
 % if n <=  10, legend('show','Location','southoutside'), end
  
  figure(2)    % plot computational full matrix eigenvalue equation relative errors 
  semilogy(errdat(:,1),errdat(:,2),'LineWidth',0.7), grid on, hold on
  a = ylim; plot([jump1,jump1],a,'k','LineWidth',1,'Displayname', ['jump1 at ',num2str(jump1)]), 
  plot([jump2,jump2],a,'k','LineWidth',1,'Displayname', ['jump2 at ',num2str(jump2)] ), hold off, 
  title('A(t) V(t) - V(t) D(t) rel. errors (via 7ifd difference equations)'),
  xlabel([' 0 \leq t \leq ',num2str(tf), ',  \eta = ',num2str(eta),',  \tau = ',num2str(tau)])
  
  figure(3)       % plot individual eigenvalue equation errors of A(t) for all t
  for jj = 1:n 
    semilogy(relerrinddat(:,1), relerrinddat(:,jj+1),'LineWidth',0.7, 'Displayname', ['eigenvalue ',num2str(jj)])
    hold on, grid on
  end  
  hold off,  
  title(['Individual equation A(t) x_k(t) - \lambda_k(t) x_k(t) errors via 7ifd difference equations']), 
  xlabel([' 0 \leq t \leq ',num2str(tf),',  \eta = ',num2str(eta),',  \tau = ',num2str(tau)])
  %if n <=  10, legend('show','Location','southoutside'), end
  
  figure(4)
  semilogy(nAdot(:,1),nAdot(:,2),'Linewidth',0.7), grid on, hold off
  title('Norm of derivative matrix for A(t)')
  xlabel([' 0 \leq t \leq ',num2str(tf),',  \eta = ',num2str(eta),',  \tau = ',num2str(tau)])
  
  figure(5)
  semilogy(VcOrthErr(:,1),VcOrthErr(:,2),'Linewidth',0.7), grid on, hold on
  a = ylim; plot([jump1,jump1],a,'k','LineWidth',1,'Displayname', ['jump1 at ',num2str(jump1)]), 
  plot([jump2,jump2],a,'k','LineWidth',1,'Displayname', ['jump2 at ',num2str(jump2)] ), hold off, 
  title('   Deviation from orthogonality of eigenvector matrix for A(t) (7ifd)')
  xlabel([' 0 \leq t \leq ',num2str(tf),',  \eta = ',num2str(eta),',  \tau = ',num2str(tau)])
  
end

function [A1,A2,A3,A4,A5,A6,Z1,Z2,Z3,Z4,Z5,Z6,jj,relerrinddat,errdat,evaluedat] = ...
              starting(kk,alt,tau,zeich,relerrinddat,errdat,evaluedat)
      
  if alt == 0, A1 = f((kk-1)*tau); else, A1 = fj((kk-1)*tau); end
  [V,D] = eig(A1); lam = diag(D)'; Z1 = [V;lam];  %sA1 = size(A1)
  if zeich == 1
    for k2 = 1:n
      relerrinddat(kk,[1,k2+1]) = [(kk-1)*tau,norm(A1*V(:,k2) - V(:,k2)*lam(k2))/abs(lam(k2))];
    end
    errdat(kk,:) = [(kk-1)*tau,norm(A1*V - V*D)/norm(A1)]; 
    evaluedat(kk,:) = [(kk-1)*tau,lam]; 
  end
  
  if alt == 0, A2 = f(kk*tau); else, A2 = fj(kk*tau); end
  [V,D] = eig(A2); lam = diag(D)'; Z2 = [V;lam];
  if zeich == 1
    for k3 = 1:n
      relerrinddat(kk+1,[1,k3+1]) = [kk*tau,norm(A2*V(:,k3) - V(:,k3)*lam(k3))/abs(lam(k3))];
    end
    errdat(kk+1,:) = [kk*tau,norm(A2*V - V*D)/norm(A2)]; % stored rowwise; [time, 1 value]
    evaluedat(kk+1,:) = [kk*tau,lam];                    % stored rowwise: [time, n evalues]
  end
  
  if alt == 0, A3 = f((kk+1)*tau); else, A3 = fj((kk+1)*tau); end
  [V,D] = eig(A3); lam = diag(D)'; Z3 = [V;lam];
  if zeich == 1
    for k4 = 1:n
      relerrinddat(kk+2,[1,k4+1]) = [(kk+1)*tau,norm(A3*V(:,k4) - V(:,k4)*lam(k4))/abs(lam(k4))];
    end
    errdat(kk+2,:) = [(kk+1)*tau,norm(A3*V - V*D)/norm(A3)];  
    evaluedat(kk+2,:) = [(kk+1)*tau,lam];
  end
  
  if alt == 0, A4 = f((kk+2)*tau); else, A4 = fj((kk+2)*tau); end
  [V,D] = eig(A4); lam = diag(D)'; Z4 = [V;lam];
  if zeich == 1
    for k5 = 1:n
      relerrinddat(kk+3,[1,k5+1]) = [(kk+2)*tau,norm(A4*V(:,k5) - V(:,k5)*lam(k5))/abs(lam(k5))];
    end
    errdat(kk+3,:) = [(kk+2)*tau,norm(A4*V - V*D)/norm(A4)];  
    evaluedat(kk+3,:) = [(kk+2)*tau,lam];
  end
       
  if alt == 0, A5 = f((kk+3)*tau); else, A5 = fj((kk+3)*tau); end
  [V,D] = eig(A5); lam = diag(D)'; Z5 = [V;lam];
  if zeich == 1
    for k6 = 1:n
      relerrinddat(kk+4,[1,k6+1]) = [(kk+3)*tau,norm(A5*V(:,k6) - V(:,k6)*lam(k6))/abs(lam(k6))];
    end
    errdat(kk+4,:) = [(kk+3)*tau,norm(A5*V - V*D)/norm(A5)]; % stored rowwise; [time, 1 value]
    evaluedat(kk+4,:) = [(kk+3)*tau,lam];                    % stored rowwise: [time, n evalues]
  end
 
  if alt == 0, A6 = f((kk+4)*tau); else, A6 = fj((kk+4)*tau); end
  [V,D] = eig(A6); lam = diag(D)'; Z6 = [V;lam];
  if zeich == 1
    for k7 = 1:n
      relerrinddat(kk+5,[1,k7+1]) = [(kk+4)*tau,norm(A6*V(:,k7) - V(:,k7)*lam(k7))/abs(lam(k7))];
    end
    errdat(kk+5,:) = [(kk+4)*tau,norm(A6*V - V*D)/norm(A6)];  
    evaluedat(kk+5,:) = [(kk+4)*tau,lam];
  end
      
%   if alt == 0, A7 = f((kk+5)*tau); else, A7 = fj((kk+5)*tau); end
%   [V,D] = eig(A7); lam = diag(D)'; Z7 = [V;lam];
%   if zeich == 1
%     for k8 = 1:n
%       relerrinddat(kk+6,[1,k8+1]) = [(kk+5)*tau,norm(A7*V(:,k8) - V(:,k8)*lam(k8))/abs(lam(k8))];
%     end
%     errdat(kk+6,:) = [(kk+5)*tau,norm(A7*V - V*D)/norm(A7)];  
%     evaluedat(kk+6,:) = [(kk+5)*tau,lam];
%   end
%      
%   if alt == 0, A8 = f((kk+6)*tau); else, A8 = fj((kk+6)*tau); end
%   [V,D] = eig(A8); lam = diag(D)'; Z8 = [V;lam];
%   if zeich == 1
%     for k9 = 1:n
%       relerrinddat(kk+7,[1,k9+1]) = [(kk+6)*tau,norm(A8*V(:,k9) - V(:,k9)*lam(k9))/abs(lam(k9))];
%     end
%     errdat(kk+7,:) = [(kk+6)*tau,norm(A8*V - V*D)/norm(A8)];  
%     evaluedat(kk+7,:) = [(kk+6)*tau,lam];       
%   end
  
  jj = kk+7; 
end

end

