a1 = 2;
a2 = 4;
a3 = 2;
b = 6;
k = 0.1;
Ti = 0.496;
Tdp = 0.1;
Tdl = 10;
L = 100;
d = (Tdl-Tdp)/L;
W2 = zeros(L,2);
Td = Tdp;
for i=1:L
 sim('lab3_P',100);
 e = simout1.data;
 t = simout1.time;
 einf = mean(e(end-10:end));
 ep1 = e - einf;
 ep2 = ep1.^2;
 tep2 = t.*ep2;
 Q = trapz(t,tep2);
 W2(i,1) = Td;
 W2(i,2) = Q;
 Td = Td + d;
end
qmin = min(W2(:,2));
Val2 = W2(:,2);
ind2 = find(Val2==qmin);
Td = W2(ind2,1);
sim('lab3_P',1000);
plot(simout1);
disp(Td);