        
function P = fun_polynomial(W, c, rankH3)

    global N
    
    w_polynomial = [];
    shift_columns = [];
    
    for idxp_w = 0:N
        w_polynomial = [w_polynomial, W^idxp_w];
        shift_columns = [shift_columns,idxp_w];
    end
    
    c = diag(w_polynomial)*c;
    c = cell2mat(arrayfun(@(x) circshift(c(x,:),[1 shift_columns(x)]),(1:numel(shift_columns))','un',0));
    
    aux_P = flip(sum(c)); % Array of the coefficients which define the polynomial in z of det[w-H]
    aux = find(aux_P~=0);
    idx_zero = aux(1);

    if rankH3 == N
        P = aux_P;
    elseif idx_zero == 1 
        idx_max_coef = N-rankH3-1;
        idx_last_coef = 2*(N+1) - idx_max_coef;

        P = aux_P(idx_max_coef:idx_last_coef);
    else 
        idx_max_coef = idx_zero;
        idx_last_coef = 2*(N+1) - idx_max_coef;

        P = aux_P(idx_max_coef:idx_last_coef);
    end


end
