function [y,x] = lsim(D,u,t,~,InterpRule)
% Linear response simulation of transfer function models.
% U is assumed to be Ns-by-Nu

%	 Author: P. Gahinet
%   Copyright 1986-2022 The MathWorks, Inc.
x = [];
if D.Ts==0 || ~(isreal(u) && isreal(D))
   % Continuous or complex 
   y = lsim(ss(D),u,t,[],InterpRule);
else
   num = D.num;
   den = D.den;
   % Computability
   D = checkCompatibility(D,'lsim',[],[],false);
   % Limit delays to simulation horizon to avoid "out of memory" errors in TFSIM
   ns = size(u,1);
   iod = min(getIODelay(D,'total'),ns);
   % Simulate
   InitState = ltipack.tfsimIC(num,den,iod);
   y = tfsim(num,den,iod,u,InitState);
end
