There is no need for the solution $\psi(x)$ to be real. What must be real is the probability density that is "carried" by $\psi(x)$. In some loose and imprecise intuitive way, you may think about a TV image carried by electromagnetic waves. The signal that travels is not itself the image, but it carries it, and you can recover the image by decoding the signal properly.
Somewhat similarly, the complex wave function that is found by solving Schrödinger equation carries the information of "where the particle is likely to be", but in an indirect manner. The information on the probability density $P(x)$ of finding the particle is recovered from $\psi(x)$ simply by multiplying it times its complex conjugate:
$$\psi(x)^*\psi(x) = P(x)$$
that gives a real function as a result. Note that it is a density: what you compute eventually is the probability of finding the particle between $x=a$ and $x=b$ as $\int_{a}^{b} P(x) dx$
As you know, when you multiply a complex number(/function) times its complex conjugate, the information on the phase is lost:
$$\rho e^{i \theta}\rho e^{-i \theta}=\rho^{2}$$
For that reason, in some places one can (not quite correctly) read that the phase has no physical meaning (see footnote), and then you may wonder "if I eventually get real numbers, why did not they invent a theory that directly handles real functions?".
The answer is that, among other reasons, complex wave functions make life interesting because, since the Schrödinger equation is linear, the superposition principle holds for its solutions. Wave functions add, and it is in that addition where the relative phases play the most important role.
The archetypical case happens in the double slit experiment. If $\psi_{1}$ and $\psi_{2}$ are the wave functions that represent the particle coming from the hole number $1$ and $2$ respectively, the final wave function is
$$\psi_{1}+\psi_{2}$$
and thus the probability density of finding the particle after it has crossed the screen with two holes is found from
$$P_{1+2}= (\psi_{1}+\psi_{2})^{*}(\psi_{1}+\psi_{2}) $$
That is, you shall first add the wave functions representing the individual holes to have the combined complex wave function, and then compute the probability density. In that addition, the phase informations carried by $\psi_{1}$ and $\psi_{2}$ play the most important role, since they give rise to interference patterns.
Comment: Feynman is quoted to have said "One of the miseries of life is that everybody names things a little bit wrong, and so it makes everything a little harder to understand in the world than it would be if it were named differently." It is quite similar here. Every book says that the phase of the wave function has no physical meaning. That is not 100% correct, as you see.
For a time-independent Hamiltonian $H$, the time-dependent Schrödinger equation $i\hbar\partial_t\Psi(x,t)=H\Psi(x,t)$ can be solved by first finding the eigenvalues und eigenstates of the Hamiltonian, that is to say by solving the time-independent Schrödinger equation $H\rho_n(x)=E_n\rho_n(x)$. Then, the time-dependent Schrödinger equation is solved by $\Psi_n(x,t)=\rho_n(x)e^{-iE_nt/\hbar}$:
$$i\hbar\partial_t\rho_n(x)e^{-iE_nt/\hbar} = E_n\rho_n(x)e^{-iE_nt/\hbar} = H\rho_n(x)e^{-iE_nt/\hbar}=H\Psi_n(x,t).$$
The general solution of the time-dependent Schrödinger equation is thus a superposition of the solutions $\Psi_n$.
Let me show you how to arrive at the above form.
Addendum 1: Separation of variables
This approach comes from the theory of differential equations, where $H$ is seen as a differential operator acting on wave functions. If the Hamiltonian is time-independent, it acts entirely only on the position variable $x$. One can then choose the separation ansatz $\Psi(x,t)=\rho(x)\phi(t)$ and substitute it into the Schrödinger equation. Dividing it by $\Psi$ gives:
$$i\hbar\frac{\partial_t\phi}{\phi} = \frac{H\rho}{\rho}$$
where on the LHS the $\rho$ could be canceled as $\partial_t$ doesn't act on it and similarly on the RHS the $\phi$ could be cancelled as $H$ only acts on the position.
Now one argues that this equation has to hold for any position $x$, so we require $H\rho/\rho$ to be constant. We call the constant $E_n$. There may be any number of solutions, thus the constant is labeled by $n$. It is of course the energy obtained from solving the time-independent Schrödinger equation: $H\rho=E_n\rho$. Replacing the RHS by the constant gives
$$i\hbar\partial_t\phi = E_n\phi\qquad\Rightarrow\qquad \phi(t)=e^{-iE_nt/\hbar}.$$
Addendum 2: Translational invariance in time The second approach is more technical but quite elegant and general. We define the time-translation operator $T_{t_0}$ whose action on a function of time is defined as $T_{t_0}f(t) = f(t-t_0)$, i.e. that function is shifted in time by $t_0$.
Clearly, the translation operator obeys a group-property $T_{t_0}T_{t_1} = T_{t_0+t_1}$. So its eigenvalues can be written in the form $f(t-t_0)=T_{t_0}f(t)=e^{\alpha t_0}f(t)$. So $f$ needs to be an exponential function, $f(t)=e^{-\alpha t}$ where $\alpha$ is arbitrary for now.
For a time-independent Hamiltonian $H$, if $\Psi_n$ is a wave function with eigenenergy $E_n$, so $H\Psi_n=E_n\Psi_n$, it holds:
$$T_{t_0}H\Psi_n(t) = T_{t_0}E_n\Psi_n(t) = E_nT_{t_0}\Psi_n(t) = E_n\Psi_n(t-t_0) = H\Psi_n(t-t_0) = HT_{t_0}\Psi_n(t),$$
so the Hamiltonian and $T_{t_0}$ commute which is denoted $[T_{t_0},H]=0$ with the commutator $[A,B]:=AB-BA$.
When two operators commute, a theorem of linear algebra states that one can find common eigenfunctions. So let's assume $\Psi_n$ is an eigenfunction to $H$ with eigenvalue $E_n$ and to $T_{t_0}$ with eigenvalue $e^{\alpha t_0}$, such that the temporal dependence is given by $\Psi\propto e^{-\alpha t}$ as reasoned before. Substituting this into the time-dependent Schrödinger equation and you get:
$$i\hbar\partial_t\Psi_n(x,t) = -i\hbar\alpha\Psi_n(x,t)\overset{!}{=}E_n\Psi_n(x,t) = H\Psi_n(x,t).$$
Therefore, $\alpha=iE_n/\hbar$, so the time-dependence of $\Psi_n$ is given by $e^{-iE_nt/\hbar}$.
This second method is kind of a Bloch theorem in time. In the usual Bloch theorem, one considers (discrete) translations in space and also finds an exponential dependence.
Best Answer
Your reasoning in the OP is correct, except when you say that the particle is at rest. Recall that in QM the exact position of a particle, and its state of motion, are not well defined. This means that it doesn't make sense to say that a particle is at rest. You could only say that a particle is at rest if $$ \Delta p\equiv 0 $$ but, because of the Heisenberg uncertainty principle, $\Delta p$ is always positive.
On the other hand, the correct statement is that it has an equal probability to be moving to the right and to be moving to the left, and thus, on average, its mean velocity is zero. You observed this in the OP, and this is the correct interpretation of $\langle p\rangle=0$.