kp=10;
kl= 0.01;
L=100;
d = (kp-kl)/L;
k = kl;


for i=0:1:L
sim('model',100)
e = simout.data;
t = simout.time;

einf = mean(e(end-20:end));
f = t.*abs(e - einf);
Q = trapz(t,f);
W(i+1,1) = k;
W(i+1,2) = Q;
k = k+d;
end

Qmin = min(W(:,2));
kindex = find(W(:,1) == Qmin);
k = W(kindex,1);

disp(Qmin);
disp(k);
sim('model',100)
plot(simout);