num = [1];
den = conv([1 3 2],[1 7 12])

T = 0.1;

sys = tf(num,den);
sysd = c2d(sys,T);

beta1 = exp(T);
beta2 = exp(2*T);
beta3 = exp(3*T);
beta4 = exp(4*T);

Amat = [1 1 1 1 1
        1 beta1 beta1^2 beta1^3 beta1^4
        1 beta2 beta2^2 beta2^3 beta2^4
        1 beta3 beta3^2 beta3^3 beta3^4
        1 beta4 beta4^2 beta4^3 beta4^4]

Bmat= [24;0;0;0;0];

DeltaU =real( inv(Amat)*Bmat)

for ind = 1:length(DeltaU)
	u(ind) = sum(DeltaU(1:ind));
end;

lsim(sysd,u)

Resp = step(sysd,3*T);
Resp1 = Resp(2);
Resp2 = Resp(3);
Resp3 = Resp(4);

e0 = 1;
e1 = 1 - DeltaU(1)*Resp1;
e2 = 1 - DeltaU(1)*Resp2- DeltaU(2)*Resp1;
e3 = 1 - DeltaU(1)*Resp3- DeltaU(2)*Resp2- DeltaU(3)*Resp1;

Gd = tf(DeltaU',conv([1 -1],[1 e1 e2 e3]),T);

tfopen = sysd*Gd
tfclose = feedback(tfopen,1);
step(tfclose)