function y = dlsym(A,B,C,D,u,t,x0) 
% y = dlsym(A,B,C,D,u,t,x0) 
% simulates the response of a discrete-time linear system to arbitrary inputs
%			
%		x(k+1) = A x(k)  +  B u(k)
%		y(k)   = C x(k)  +  D u(k)
%
%	A    :   n by n  dynamics    matrix
%	B    :   n by m  input       matrix 
%	C    :   l by n  output      matrix
%	D    :   l by m  feedthrough matrix	
%	u    :   l by p  matrix of sampled inputs
%	t    :   1 by p  vector of uniformly spaced points in time
%	x0   :   n by 1  vector of initial states, defaults to zero
%	y    :   m by p  matrix of the system outputs
%

 [n,m,l] = abcddim(A,B,C,D);

 points = size(u,2);		 % number of data points 

 if ( nargin < 6 )
   t = [1:points];
 end

 if ( nargin == 7 )
   x = x0;                       % initial conditions for the state 
 else
    x = zeros(n,1);              % initial conditions are zero
 end

 y = NaN(l,points);            % memory allocation for the output

 y(:,1) = C * x  +  D * u(:,1);

 for p = 2:points

     y(:,p) = C * x  +  D * u(:,p);

     x = A * x  +  B * u(:,p);

 end

end % --------------------------------------------------------- dlsym.m
% 2021-07-19 ...
%   replaced ...  x = A * x  +  B * u(:,p); 
%   ... with ...  x = A * x  +  B * u(:,p-1);
% 2023-10-01
%  ... switch order of dynamics and output eqn calc's


