function y = simul(G,u,err)
%simul Simulates discrete system G with input u. Works for improper
%transfer functions

num = G.Numerator{1};
den = G.Denominator{1};

p = sum(num~=0) - sum(den~=0);

t = linspace(0,length(u)-1,length(u));

if p > 0
den2 = [den zeros(1,p)];
G2 = tf(num,den2,G.Ts);
y_shift = lsim(G2,u,t);
y = [y_shift(p+1:end) ; y_shift(end)*ones(p,1)];
else
y = lsim(G,u,t);
end
if nargin == 3
y = y + err;
end
end