To talk about "the" wavefunction doesn't really make sense. There are many wavefunctions in the Hilbert space of this system for example. Usually what is meant with this language are the eigenstates of some model Hamiltonian H.
In the nearly free electron model, we consider non interacting electrons with the effective Hamiltonian
$$ H = \frac{p^2}{2m} + V $$
The potential energy comes from a combination of electron-electron and electron-nucleon interaction and can be modeled as a small perturbation. Note that we are considering a single electron problem now. Further, because of the crystal symmetry, the potential has lattice periodicity
$$V(R+a) = V(R)$$
for lattice vectors $a$. We want to find the energy eigenstates of this system. Label the electron plane waves as $\psi_k = e^{ikx}$. It is clear that those are just eigenstates of the kinetic energy operator, and not of the full Hamiltionian due to the V. The important point is now that due to the lattice periodicity, the potential energy only couples plane waves if the difference of the wavenumbes is a multiple of a reciprocal lattice vector $G=\frac{2 \pi}{a}n$. In matrix notation:
$$\langle\psi_k,V\psi_{k'}\rangle=\sum_G\delta_{k-k',G} V_G$$
(You should verify this as an exercise)
The $V_G$ falls off really fast for increasing $G$ so usually it is enough to consider the coupling of only two plane wave states. Now, if we have a small perturbation, its effect will be strongest on states that were degenerate without the perturbation. This means that we can ignore the effect of $V$ on plane waves where there is no other wave with the same energy. Therefore, the states with small $k$ inside your band are still plane waves to a good approximation. However, if you go to the edge of the band, the state with $k=\frac{\pi}{a}$ has the same energy as the one with $k=-\frac{\pi}{a}$, and they are coupled because they differ by a reciprical lattice vector. In this two state subspace the Hamiltonian looks like this:
$$H= \begin{pmatrix}
\epsilon & V_{\frac{2 \pi}{a}}\\
V_{\frac{2 \pi}{a}} & \epsilon
\end{pmatrix}$$
Because the $\epsilon$ are equal, the eigenstates are $\phi_{\pm}=\psi_{\frac{\pi}{a}} \pm \psi_{-\frac{\pi}{a}} $. (If the electron is in this state, it is in a superposition of traveling into two different directions.) Edit: Actually I think we should be careful about the velocity of the electron. In the tight binding model for example, an eigenstate of $k$ is also an eigenstate of the velocity. Note that the states $\psi_k = e^{ikx}$ are coarse grained and cannot be the eigenfunctions if we look on length scales of interatomic distances. Very close to the nucleus, V will not be a perturbation anymore. Thus the tight binding approach should yield the correct velocity.
In fact, not only the two electron states which sit exactly at the band gap are effected by the potential. Also those which lie close to the gap. The eigenstates of $H$ are still eigenfunctions of the crystal momentum $k$, but not longer eigenfunctions of the kinetic momentum $p=-i \nabla$.
In kittel's book on solid state physics it says that wavefunctions at $k=\pm\frac{\pi}a$ are not the traveling waves exp($i\pi x/a)$ or exp($-i\pi x/a)$ of free electrons . At these special values of k the wavefunctions are made of equal parts of waves traveling to the right and to the left. When the bragg reflection conditon $k=\pm\frac{\pi}a$ is satisfied by the wavevector a wave traveling to the right is Bragg-reflected to travel to the left and vice versa.
