Even with single-particle systems, analytical solutions are only possible for a small subset of quantum mechanical problems: https://en.wikipedia.org/wiki/List_of_quantum-mechanical_systems_with_analytical_solutions For example, if you were to change the external potential of the particle in a box into something even a tad weirder, like a polynomial function, I would surmise that it is impossible to solve it analytically. The same goes with the one-electron operators in Hartree–Fock theory, because you have the mean-field potential to contend with.
Therefore, we need to bring in numerical methods. One way of doing that is to express the eigenstates as a linear combination of certain basis functions. If the set of basis functions – the basis set – is complete, then in principle we can express any function we like as a linear combination of these basis functions. Mathematically, this is expressed by a "completeness relation":
$$\hat{1} = \sum_i |i\rangle\langle i| \quad \Longleftrightarrow \quad |f\rangle = \hat{1}|f\rangle = \sum_i |i\rangle\langle i | f \rangle = \sum_i c_i|i\rangle \text{ where }c_i = \langle i | f \rangle$$
Unfortunately, in the case of an atom or molecule, we need an infinite number of basis functions to have a complete set. This is obviously not possible, so we need a finite number of basis functions, judiciously chosen such that the result obtained using this limited basis is close enough to the result obtained with the infinite basis. The AOs (or approximations to them) are extremely convenient in this regard, and thus LCAO is born.
