kp = 10;
kl = 0.01;
L = 1000;
d = (kp - kl) / L;
W = zeros(L, 2);
k = kl;

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

    einf = mean(e(end-20:end));
    f = t.*(e - einf).^2;
    Q = trapz(t, f);

    W(i+1, 1) = k;
    W(i+1, 2) = Q;
    k = k + d;
end

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

disp(Qmin);
disp(k);
sim('kolos1_test');
plot(simout);