A={{1,2,4,67},{2,81,42,1},{4,42,6,-2/3},{67,1,-2/3,10}};
A//MatrixForm

(*A is symmetric*)
Transpose[A]==A

(*Characteristic polynomial *)
Simplify[Det[A-λ IdentityMatrix[4]]]

(*Coeff of order 0 and 3*)
(-1)^(4-4)Det[A]
(-1)^(4-3)Tr[A]

(*Find roots*)
NSolve[(51077416/9)+((3642910 λ)/9)-((43393 (λ^2))/9)-98 (λ^3)+(λ^4)==0,λ]

(*Find eigenvalues and corresponding eigenvectors*)
eigensys=Eigensystem[A];
(*Eigenvalues*)
{l1,l2,l3,l4}=N[eigensys[[1]]]

(*Eigenvectors*)
v1=N[eigensys[[2,1]]]
v2=N[eigensys[[2,2]]]
v3=N[eigensys[[2,3]]]
v4=N[eigensys[[2,4]]]

(*Checking that these are indeed eigenvectors*)
A.v1
l1 v1

A.v2
l2 v2

(*Normalizing eigenvectors to be unit length*)
u1=v1/Sqrt[v1.v1];
u2=v2/Sqrt[v2.v2];
u3=v3/Sqrt[v3.v3];
u4=v4/Sqrt[v4.v4];
P=Transpose[{u1,u2,u3,u4}];
P//MatrixForm

Inverse[P]//MatrixForm

Transpose[P]//MatrixForm

(*So P is orthogonal! Indeed, symmetric matrices are orthogonally diagonalizable.*)

(*Watch out for small numerical errors*)
N[Transpose[P]-Inverse[P]]
Inverse[P]==Transpose[P](*This might evaluate to False, depending on your setup*)