Here is Magma code that will construct a 10x10 example, with one eigenvalue about +10 and nine of -1. Note that Magma has little numerical linear algebra, so I do it myself.
As can be seen in the final $P$ and $Q$, the numbers become of size about $10^{40}$ here. You can round to the nearest integer after multiplying by $10^{100}$ and get an integral answer of course. It seems much harder to require that $P$ and $Q$ have positive coefficients.
RF:=RealField(100);
function ThetaRoot(M)
h11,h12,h21,h22:=Explode(Eltseq(M));
t2:=h11^2; t1:=h22^2; t12:=(h12+h21)^2-2*h11*h22;
a:=t1+t2+t12; b:=-2*t1-t12; c:=t1;
R:=Roots(Polynomial([c,b,a]),RF); if #R eq 0 then R:=[<-b/2/a,2>]; end if;
th:=Arccos(Sqrt(R[1][1])); return th; end function;
function UnitarilySimilar(M)
t:=Trace(M); d:=Degree(Parent(M)); S:=SymmetricGroup(d);
M:=M-t/d*Parent(M)!1;
if d eq 1 or M[1][1] eq 0 then return Parent(M)!1; end if;
D:=Diagonal(M);
POS:=[i : i in [1..d] | D[i] gt 0]; NEG:=[i : i in [1..d] | D[i] lt 0];
p:=POS[1]; n:=NEG[1]; u:=S!1;
if p ne 1 then u*:=S!(p,1); end if; if n ne 2 then u*:=S!(n,2); end if;
if p eq 2 and n eq 1 then u:=S!(1,2); end if;
T1:=PermutationMatrix(BaseRing(M),u);
M1:=T1*M*Transpose(T1);
if Abs(M1[1][1]) gt Abs(M1[2][2]) then
Tt:=Parent(M)!1; Tt[1][1]:=0; Tt[1][2]:=-1; Tt[2][1]:=1; Tt[2][2]:=0;
T1:=Tt*T1; M1:=T1*M*Transpose(T1); end if;
th:=ThetaRoot(Submatrix(M1,[1,2],[1,2]));
T2:=Parent(M)!1;
T2[p][p]:=Cos(th); T2[p][n]:=-Sin(th);
T2[n][p]:=-T2[p][n]; T2[n][n]:=T2[p][p];
M2:=T2*M1*Transpose(T2);
if Abs(M2[1][1]) gt 10^(-25) then T2[p][n]:=Sin(th); T2[n][p]:=-T2[p][n];
M2:=T2*M1*Transpose(T2); end if;
U:=UnitarilySimilar(Submatrix(M2,2,2,d-1,d-1));
U:=DirectSum(<DiagonalMatrix([BaseRing(U)!1]),U>);
return U*T2*T1; end function;
M:=DiagonalMatrix([RF!10,-1,-1,-1,-1,-1,-1,-1,-1,-1]); d:=Degree(Parent(M));
U:=UnitarilySimilar(M); MT:=U*M*Transpose(U); I:=Parent(MT)!1;
for a in [1..d] do for b in [a+1..d] do // Ballantine gives "too diagonal"
I[a][b]:=Random([-2^25..2^25])/2^32; end for; end for; // perturb it
PERTURB:=I*MT*Transpose(I); B:=PERTURB; TRANS:=I*U;
Diagonal(PERTURB);
for a in [1..d] do for b in [a+1..d] do B[a][b]:=0; end for; end for;
for a in [1..d] do B[a][a]:=B[a][a]/2; end for;
function EpsilonKernel(M) d:=Degree(Parent(M));
S:=Vector([BaseRing(M)!0 : i in [1..d]]);
for e in [1..d] do T:=M;
for f in Reverse([1..e-1]) do S[f]:=T[e][f]/T[f][f];
T[e]-:=T[f]*S[f]; end for;
if Abs(T[e][e]) lt 10^(-50) then S[e]:=-1; return S/Sqrt(Norm(S)); end if;
end for; end function;
function MagmaBrainDeadEigenvectors(M)
R:=[r[1] : r in Roots(CharacteristicPolynomial(M))];
return DiagonalMatrix(R),Matrix([EpsilonKernel(M-r*Parent(M)!1) : r in R]);
end function;
D,V:=MagmaBrainDeadEigenvectors(B);
H:=Transpose(V)*D*V; P:=Transpose(B)*H; Q:=H^(-1);
P:=(P+Transpose(P))/2; Q:=(Q+Transpose(Q))/2; // cheat for symmetry
Roots(CharacteristicPolynomial(P*Q+Q*P));
EDIT: In larger dimensions, the $Q=H^{-1}$ becomes time-consuming, perhaps for stability.