function passed = check_ILC_statetransitionmatrix(verbosity)

if nargin==0
  verbosity=1;
end

tol = 1e-9; % numerical tolerance

nx = 5;
nt = 21;
ny = 2;
nu = 2;

for ik = 1:nt
  b=randn(nx);
  A{ik} = 0.99*(b*b')/max(eig(b*b')); % hermitian with eigenv <0.99
  B{ik} = ones(nx,nu);
  C{ik} = randn(ny,nx);
  D{ik} = zeros(ny,nu);
end
sys.A = A; sys.B=B; sys.C=C; sys.D=D;

iT = 1:nt;
U = randn(nu,nt);
x0 = randn(nx,1);

%
[Hk,Tk] = ILC_statetransition_matrix(sys);

%%
dU = randn(size(U));
dx0 = 1e-5*rand(nx,1);

dY1 = ltvsim(sys,U+dU,x0+dx0)-ltvsim(sys,U,x0);
%
du = reshape(dU,1,[])';
dy2 = Tk*du + Hk*dx0;
dY2 = reshape(dy2',ny,[]);

pass(1) = (norm(dY2-dY1)/norm(dY1) < tol);

if verbosity>0
  figure(9)
  subplot(221)
  plot(iT,dY1,'b',iT,dY2,'r--')
  title('dY')
  ylabel('dY')
  
  subplot(222)
  plot(iT,dY2-dY1);
  title('error')
end

%% Check reduced time grids
ity = round(nt/3):nt;
itu = 1:round(2*nt/3);

% get transition matrices only for cetain time points
[Hk,Tk] = ILC_statetransition_matrix(sys,ity,itu);

% perturb only those inputs
dU2 = zeros(size(dU)); dU2(:,itu) = dU(:,itu); 
dYY = ltvsim(sys,U+dU2,x0+dx0)-ltvsim(sys,U,x0);

dY1 = dYY(:,ity);

du = reshape(dU(:,itu),1,[])';
dy2 = Tk*du + Hk*dx0;
dY2 = reshape(dy2',ny,[]);

pass(2) = (norm(dY2-dY1)/norm(dY1) < tol);

if verbosity>0
  figure(9)
  subplot(223)
  plot(iT(ity),dY2,'b',iT(ity),dY1,'r--')
  title('dY')
  ylabel('dY on reduced time grid')
  subplot(224)
  plot(iT(ity),dY2-dY1);
  title('error');
  
end

%%


%%
passed = all(pass);
if verbosity>0
  if passed
    disp('passed');
  else
    disp('failed');
  end
end