I solved equation $H[v] = eS[v]$ using some numerical package (e,V = scipy.lingalg.eig(H,S))
There $[v]$ is matrix of all generalized eigenvectors $v_i$, $H$ is projection of operator $\hat H$ on non-orthogonal normalized basis functions $[\chi]$
$H_{ij} = \int \chi_i \hat H \chi_j$
and $S = [\chi]^T[\chi]$ is overlap matrix
$S_{ij} = \int \chi_i \chi_j$
The matrix of generalized eigenvectors $[v]$ is normalized, but it is not ortho-normal
According to wikipedia $[v]$ it should be S-orthogonal, meaning $[v]^TS[v]$ is diagonal, but not ortho-normal. (that is right, checked it).
But what I want and what I mean by true eigenvector is unitary matrix $U = [u]$ such that it transforms the non-orthogonal basis functions $[\chi]$ to orthonormal functions $[\phi]$ which diagonalize operator $\hat H$
$[\phi] = U^T [\chi] $
$\hat H \phi_i = \lambda_i \phi_i $
NOTE:
$U$ can be obtained by Löwdin tranform, that is solvin $(S^{-1/2} H S^{-1/2})u=\lambda u$ rather than $Hv=eSv$, but I don't want to do that since to obtain $S^{-1/2}$ I need to diagonalize $S$-matrix, that is do two diagonalizations rather than one.