Chemistry - Linear Combination of Atomic Orbitals
Solution 1:
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 "projectively". 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: $$ <\Psi_i\vert H\vert\Psi_j> = \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$: $<\Psi\vert H\vert\Psi>$, the "expectation value" 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.
Solution 2:
As you requested a technically simple answer, I'll describe every step thoroughly to be safe. I'm also assuming we're dealing with the time-independent Schrödinger equation, and will refer to it below just as the Schrödinger equation (I'm happy to explain why in the comments if you're interested).
By solving the Schrödinger equation for your system, you find the eigenstates/eigenfunctions of your hamiltonian, as well as the energy values of each of those eigenstates. Because your state function (the specific wave function you're interested in -- the one that describes your system in its current state) obeys the Hamiltonian of the system, then you are guaranteed that the state function can be exactly written (with no loss of accuracy or approximation) as a linear combination of the eigenfunctions you found by solving the Schrödinger equation.
For any system containing two or more electrons, the electron-electron interaction terms render it impossible to separate the Schrödinger equation into separable variables, meaning you'll need to use further approximations to solve the problem. However, within the Born-Oppenheimer approximation, the Schrödinger equation can be solved exactly when there's a single electron. That being said, the dihydrogen cation is a great example of where that doesn't mean it can be done easily.
A linear combination of atomic orbitals is a workaround to having to solve the Schrödinger equation for your system. Using them is similar in mindset to using the eigenfunctions to compose your state function. However, this time, there is no guarantee that your wave function of interest will exactly be a mixture of the eigenfunctions. This is because here, your eigenfunctions are for the hydrogen atoms, which is not exactly the system we are looking at! As such, you often need to use numerous basis functions to get closer to agreement with the exact solution. The usage of this substitution will certainly come with an optimization through a self consistent field or other theoretical method, in order to find the best possible combination of the eigenfunctions to reflect your state function (you wouldn't expect a 2p and a 1s atomic basis function to equally contribute to the dihydrogen ground state wave function).
TLDR: The atomic orbitals are not eigenfunctions of the dihydrogen cation, and so substituting the actual state function with a linear combination of atomic orbitals will be approximate. However, if you use the complete basis set (CBS) limit and use a sufficient level of theory, you can recover excellent agreement with the exact solution.
Solution 3:
I am paraphrasing from an article in Wikipedia and applying it to your specific case.
To simplify things and to obtain a method more readily generalizable to more complex species, we use the linear combination of atomic orbitals concept.
The special thing about the ground state solution of the Schrödinger equation is that it has the lowest energy. The variational method uses the idea that any other wave function will have a higher energy.
If you have a bunch of wave functions and want to know which one is closest to the real one for the ground state, you calculate the energy each one of these would have (using the Hamiltonian for $\ce{H2+}$ in your case). Instead of taking unrelated functions, you can also build trial wave functions from linear combinations of basis functions, in your case just the 1s functions of the two hydrogen atoms, each multiplied by a coefficient. This is called a linear combination. Then, you optimize those two coefficients until you get the lowest energy.
How do we know that the wave function obtained from the linear combination obtained (say of two 1s orbitals of each atom) satisfies the Hamiltonian of this new system involving two nuclei?
They don't satisfy the Hamiltonian exactly (i.e. they are not eigenfunctions of the Hamiltonian) but they come pretty close. We know that by comparing the calculated energy to experimental data or to more sophisticated computational results. In the case of the $\ce{H2+}$ ion, you can also compare it to the exact solution.
With the 1s function centered on each atom as basis functions, you get a pretty good solution (with two positive coefficients) and a pretty awful solution (with coefficients of opposite signs), the latter of which is called the anti-bonding case (i.e. an excited state that does not hold the atoms together).
Do we need to make some iterative calculation as in the Hartree-Fock self consistent field method.
No, you don't need to iterate because there is only one electron, so there is no need for a potential energy coming from the other electrons (i.e. the "self-consistent field"). Also, with only two coefficients to vary, you will be able to find the global minimum in a straightforward manner.
If you want better solutions, you add more basis functions to your trial function (like 2s, 2p etc, or something completely different). The calculation will get more complicated, but the solutions will be closer to the actual ground state wave function. Now it will be more difficult to find the set of coefficients that give you the wave function with the lowest energy, and there might be some iteration to search for the global minimum, but that iteration is different from the ever-changing Hamiltonian in the Hartree-Fock method.
There are some tricks to cut down on computational time (or to do some of this qualitatively in your head or on paper). One central idea is to sort basis functions by symmetry, and keep different symmetries separate from each other (certain terms will be zero when the two basis functions appearing in the are zero). Another trick is to pre-compute parts of some terms, and then just vary the coefficients.
In your case no tricks are necessary because you are just combining two basis functions.