function [ K,S ] = solveGains( lambda, LB, A, C, Q, R )

[outs,states] = size(C);
S = zeros(states,states,LB);                     % matrici per il calcolo ricorsivo
K = zeros(states,outs,LB);                        % guadagni statici (N)

lambda = [lambda lambda(end)*ones(1,LB-length(lambda))];

% soluz asintotica riccati modificata
for j = 1:100
    S(:,:,LB) = A*S(:,:,LB)*A' + Q - lambda(LB)*A*S(:,:,LB)*C'*pinv(C*S(:,:,LB)*C' + R)*C*S(:,:,LB)*A';
end

% guadagno ultima posizione buffer
K(:,:,LB) = S(:,:,LB)*C'*pinv(C*S(:,:,LB)*C' + R);

% a scendere calcola il guadagno fino a posizione buffer 1
for j = (LB-1):-1:1
    S(:,:,j) = lambda(j)*A*(eye(4)-K(:,:,j+1)*C)*S(:,:,j+1)*(eye(4)-K(:,:,j+1)*C)'*A' +...
        (1-lambda(j))*A*S(:,:,j+1)*A' + Q + A*K(:,:,j+1)*R*K(:,:,j+1)'*A'; 
    K(:,:,j) = S(:,:,j)*C'*pinv(C*S(:,:,j)*C' + R);
end
