function [Y,X] = ltvsim(sys,U,x0)

assert(iscell(sys.A),'A not a cell');
assert(iscell(sys.B),'B not a cell');
assert(iscell(sys.C),'C not a cell');
assert(iscell(sys.D),'D not a cell');

nx = numel(x0);
ny = size(sys.C{1},1);

nt = numel(sys.A);

X = zeros(nx,nt);
Y = zeros(ny,nt);

X(:,1) = sys.A{1}*x0 + sys.B{1}*U(:,1);
Y(:,1) = sys.C{1}*X(:,1)   + sys.D{1}*U(:,1);

for ik=2:nt;
    X(:,ik)   = sys.A{ik}*X(:,ik-1) + sys.B{ik}*U(:,ik);
    Y(:,ik)   = sys.C{ik}*X(:,ik)   + sys.D{ik}*U(:,ik);
end

return