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.