Mixing of molecular orbitals in polyatomic molecules I am trying to solve the following Hamiltonian,$$H=-\frac{1}{2}\nabla^{2}-\sum_{i=1}^{10}\frac{1}{|\vec{r}-\vec{R_{i}}|}$$ Here $Z=1$ . To get all ther molecular orbitals of this system. I took the following linear combination to construct the hamiltonian matrix,$$\psi(\vec{r})=\sum_{i=1}^{10}c_{i}1s(\vec{r}-\vec{R_{i}})+\sum_{i=1}^{10}a_{i}2p_{x}(\vec{r}-\vec{R_{i}})+\sum_{i=1}^{10}b_{i}2p_{y}(\vec{r}-\vec{R_{i}})+\sum_{i=1}^{10}g_{i}2p_{z}(\vec{r}-\vec{R_{i}})$$

I took this wave functions and solved for $E$ using the following equation. Usually in polyatomic molecules, one would look for symmetry to construct various molecular atomic orbitals. However, here there is no symmetry. How do I get the molecular orbitals? What combination should I have taken to solve this problem. Please suggest some references where I can read the theory of molecular orbitals for polyatomic molecules with no symmetry
The standard way is to solve for the eigenvectors of your Fock matrix. Since you have a single electron, that's your Hamiltonian. Any number of programming languages have algorithms for this implemented. Each of your coefficients would correspond to an entry in the eigenvector.