Since one expects the ground state to share the same symmetry as the potential, a reasonable guess for the ground state would be a gaussian; from there you could try f(x)exp[-x^2/2]. That's certainly a valid method of solving the equation. However, as you're looking for something more mechanical (the following comes almost verbatim from Mike Stone/Paul Goldbart's notes for the Math Methods course at UIUC):
Consider the operator \hat{\mathcal{H}} = -\partial_x^2 + x^2. The Harmonic Oscillator ODE is then \hat{\mathcal{H}}\psi = E\psi.
Consider also the operators Q = \partial_x +x and Q^\dagger = -\partial_x + x and notice that \mathcal{H} = Q^\dagger Q + 1 and \mathcal{H} = QQ^\dagger - 1.
Suppose \phi is an eigenfunction of Q^\dagger Q with eigenvalue \lambda. It follows that Q\phi is an eigenfunction of QQ^\dagger with the same eigenvalue:
Q^\dagger Q\phi = \lambda \phi
Now apply the operator Q to both sides:
QQ^\dagger (Q\phi) = \lambda (Q\phi).
Now, there are two ways Q\phi could fail to be an eigenfunction. One is that it is the zero function, but then this means that the LHS is zero and hence the eigenvalue was zero too. Conversely, the eigenvalue could be have been zero to start with. But then this means that the inner product \left<\phi,Q^\dagger Q\phi\right> = \left<Q\phi,Q\phi\right> = 0, and hence Q\phi was zero. Accordingly, QQ^\dagger and Q^\dagger Q have the same spectrum except for any possible zero eigenvalues. Now, of course, there are zero eigenvalues. Solving the zero eigenvalue problems, you see that Q^\dagger Q \phi = 0 is solved by \phi_0 = e^{-xx^2/2}, which is normalizable. The other ordering of the pair also has a zero eigenvalue solution, but it is not normalizable and so you throw it out.
Now, using the relation between H and the Q's, you see that \phi_0 is an eigenfunction of H with eigenvalue 1, and an eigenfunction of QQ^\dagger with eigenvalue 2. Accordingly, Q^\dagger \phi_0 is an eigenfunction of Q^\dagger Q with eigenval 2 and H with eigenval 3. Keep iterating this process to get \phi_n = (Q^\dagger)^n \phi_0 is an eigenfunction of H with eigenvalue 2n + 1. I guess it turns out you can write Q^\dagger = -e^{x^2/2}\partial_x e^{-x^2/2}, and in that way you generate the Hermite Polynomials.
So that's one way to determine the solution 'mechanically'.