Linear Combination of Atomic Orbitals
Using a Linear Combination of Atomic Orbitals (LCAO) will not exactly satisfy the Schrödinger Equation / be an eigenfunction of the Hamiltonian
We can write down the Hamiltonian for the $\ce{H2+}$ cation with bond length $2R$, nuclei at $\vec{R_1}$ and $\vec{R_2}$, using the Born-Oppenheimer approximation:
$$
H\Psi(\vec{r}) = \left[-\frac{\nabla^2 }{2} - \frac{1}{\vert\vec{r}-\vec{R_1}\vert} - \frac{1}{\vert\vec{r}-\vec{R_2}\vert} \right]\Psi(\vec{r}) = E \Psi(\vec{r})
$$
If you feed in a linear combination of 2 $1s$ orbitals into $H\Psi$ you get a mess that is not a multiple of the linear combination you put in - i.e. a linear combination of 1s orbitals is not the exact solution.
In this case, because we only have one electron, the Schrödinger Equation can be solved to arbitrary accuracy using a variety of numerical techniques for second order partial differential equations in 3D; such as solution in an basis of functions or on a grid in space.
Solving in Practice
A suitable infinite basis of functions can describe all possible functions in 3D.  Two different examples are spherical harmonics with suitable radial components (Slater type orbitals) such as the $1s$ orbital, or Gaussian functions with suitable angular components (Gaussian type orbitals more commonly used in computational chemistry today).
With a infinite number of functions we would have an exact solution, but with only the first few terms in a suitable series we can be remarkably accurate.
In practical electronic structure theory, this means 3-30 functions per atom with suitable chosen parameters, centering the functions on atomic nuclei, as that's where the electrons will be mostly due to attraction to nuclei.
Variational Principle
Seeing as using a finite number of orbitals will not give an eigenfunction of the Hamiltonian, and certainly not the lowest energy ground state, we need another way to confirm a solution.
Instead we solve the Schrödinger equation &quot;projectively&quot;.  Hamiltonians are hermitian operators.  One of the properties of hermitian operators is that their eigenfunctions are  orthogonal.  What this means is that for a Hamiltonian $H$ and two exact solutions $\Psi_i$ and $\Psi_j$ the integral:
$$
&lt;\Psi_i\vert H\vert\Psi_j&gt; = \iiint_{\text{all space}} \Psi_i^* H \Psi_j \;\text{d}V = \cases{E_i \text{ if } \Psi_i = \Psi_j\\0 \text{ otherwise}} 
$$
Taking the integral of two wavefunctions, one multiplied by the complex conjugate of the other, is often referred to as projecting one wavefunction onto another.
Another property of the Hamiltonian operator is that the value of the integral with some generalised wavefunction $\Psi$:
$&lt;\Psi\vert H\vert\Psi&gt;$, the &quot;expectation value&quot; of the energy of $\Psi$, is minimised by the exact ground state wavefunction.  This property is know as the variational principle.  Hence in self consistent field/Hartree Fock calculations, we solve equations to minimise that integral of our LCAO wavefunctions by varying the coefficents to get the best LCAO approximation to the true wavefunction.