The meaning of the phase in the wave function

This is an important question. You are correct that the energy expectation values do not depend on this phase. However, consider the spatial probability density $|\psi|^{2}$. If we have an arbitrary superposition of states $\psi = c_{1} \phi_{1} + c_{2} \phi_{2}$, then this becomes

$|\psi|^{2} = |c_{1}|^{2}|\phi_{1}^{2} + |c_{2}|^{2} |\phi_{2}|^{2} + (c_{1}^{*} c_{2} \phi_{1}^{*} \phi_{2} + c.c.)$.

The first two terms do not depend on the phase, but the last term does. ($c_{1}^{*}c_{2} = |c_{1}||c_{2}|e^{i (\theta_{2} - \theta_{1})}$). Therefore, the spatial probability density can be heavily dependent on this phase. Remember, also, that the coefficients (or the wavefunctions, depending on which "picture" you are using) have a rotating phase angle if $\phi_{1,2}$ are energy eigenstates. This causes the phase difference $\theta_{2} - \theta_{1}$ to actually rotate at the energy difference, so that $|\psi|^{2}$ will exhibit oscillatory motion at the frequency $\omega = (E_{2} - E_{1})/\hbar$. This is known as Rabi oscillation, and is also related to optical transitions and many other quantum phenomena.

In summary, the phase information in a wavefunction holds information, including, but not limited to, the probability density. In a measurement of energy this is not important, but in other measurements it certainly can be.

Also you can modify the wavefunction with a global phase $\psi(x)\rightarrow e^{i\phi}\psi(x)$ without affecting any expectation values because the phase factor will cancel when taking inner products, so this global phase doesn't contain any information. Only relative phases are meaningful in quantum mechanics.

For a particle of mass $m$ with a simple Hamiltonian in position-space $\mathcal{H} = -\frac{\hbar^2}{2m}\nabla^2+V(\vec{x})$, if you write a general wavefunction as $$\Psi(t;\vec{x}) = \sqrt{\rho}e^{iS/\hbar}\text{,}$$ where $S$ and $\rho\geq 0$ are real, then the phase information $S$ directly corresponds to the probability current $$\mathbf{J} = \frac{\rho}{m}\nabla S\text{,}$$ the continuity equation for which turns out to be exactly the imaginary component of the Schrödinger equation, $$\frac{\partial\rho}{\partial t} + \nabla\cdot\mathbf{J} = 0\text{.}$$ As expected from more general considerations, a global phase factor is irrelevant because only its gradient appears. As a side note, the real component of the Schrödinger equation turns out to be the classical Hamilton–Jacobi equation corrected by one extra term proportional to $\hbar^2$.

The probability current can also be defined in more complicated situations, but it remains the case that morally speaking, the phase information is critical to how the wavefunction evolves in time.