function [passed, norm_diff] = check_ILC_Tk(ILC,out0,model,params,x0,g0,v0,U,outfun,varargin)

simres_arg = false;
if numel(varargin)>0
  for ii = 1:numel(varargin)
    if isfield(varargin{ii},'Xdot')
      simres0 = varargin{ii};
      simres_arg = true;
    else
      verbosity = varargin{ii};
    end
  end
end

%% optional check of Tk
dU = zeros(size(U));
nu = numel(ILC.iu);
nt = size(U,2);

% generate dU perturbation
pert = ones(nu,1) * 5*sin((0:nt-1)/nt*5);
dU(ILC.iu,ILC.itu) = pert(:,ILC.itu);
inegU = (U+dU)<0 ; % indices of negative U
dU(inegU) = -dU(inegU); % flip

% find corresponding p
du = reshape(dU(ILC.iu,ILC.itu),[],1);
dp = (ILC.dudp\du);

% recalculate U
dup = ILC.dudp*dp;
dU(ILC.iu,ILC.itu) = reshape(dup,nu,[]);

% perturbed sim
simres2 = RAPTOR_predictive(x0,g0,v0,U+dU,model,params);
out2 = RAPTOR_out(simres2,model,params);

if simres_arg
  dY2 = outfun(simres2)-outfun(simres0);
else
  dY2 = outfun(out2)-outfun(out0);
end
if ~isempty(ILC.iy)
  dY2 = dY2(ILC.iy);
end
% linearization
dy = ILC.TTk*dp;
dY = reshape(dy,ILC.ny,[]);

%% norms
norm_diff = (norm(dY2-dY)/norm(dY));
tol = 1e-4;
passed = (norm_diff < tol);

if verbosity>0
  if ~passed
  fprintf('\nILC Tk matrix test failed because too big difference between linearized and real model %.3e>%.3e.\n',norm_diff,tol)
  else
    fprintf('\nILC Tk matrix test passed\n');
  end
  
%% optional plots
if verbosity>0
  %%
  ttu = out0.time(ILC.itu);
  tty = out0.time(ILC.ity);
  
  figure(99); clf;
  subplot(311)
  plot(ttu,U(ILC.iu,ILC.itu),'b'); hold on;
  plot(ttu,U(ILC.iu,ILC.itu)+dU(ILC.iu,ILC.itu),'r--')
  legend('original trajectory', 'perturbed trajectory')
  xlabel('time')
  
  if numel(tty)>1
    % time plot
    subplot(312)
    hp(1,:)=plot(tty,dY2,'b'); hold on;
    hp(2,:)=plot(tty,dY,'r--');
    legend(hp(:,1),'full model', 'linearised model')
    title('perturbation')
    xlabel('time')
    subplot(313)
    plot(tty,dY2-dY,'b');
    title('error');
    xlabel('time')
  else
    % profile plot
    subplot(312)
    iiy = 1:numel(dY2);
    hp(1,:)=plot(iiy,dY2,'b'); hold on;
    hp(2,:)=plot(iiy,dY,'r--');
    legend(hp(:,1),'full model', 'linearised model')
    title('perturbation')
    xlabel('index')
    subplot(313)
    plot(iiy,dY2-dY,'b');
    title('error');
    xlabel('index');
  end
  
end
end
%%
