The Newton's equation of motion in potential $V(x)$ are:
$$
\dot{x}=\frac{p}{m},\dot{p}=-\frac{dV(x)}{dx}
$$
or simply
$$m\ddot{x}=-\frac{dV(x)}{dx}=F(x),$$
where $F(x)= -\frac{dV(x)}{dx}$ is called (conservative) force.
If we deal with a quantum mechanical particle, its momentum and position are related via Heisenberg uncertainty relation, which means that the joint equations of $x,p$, known simultaneously to arbitrary precision, such as the Newton equations, are meaningless. However, the correspondence principle suggests that, if we were dealing with a narrow wave packet, we would still expect write the Newton equations for the averages of position and momentum, so that classical mechanics emerges smoothly from the quantum one:
$$
\frac{d}{dt}\langle x\rangle=\frac{\langle p\rangle}{m},
\frac{d}{dt}\langle x\rangle=-\frac{dV(\langle x\rangle)}{dx}
$$
Ehrenfest theorem is the first step in this direction. Indeed, for Hamiltonian $$H=\frac{p^2}{2m}+V(x)$$ we could write Heisenberg equations of motion (formally identical to Newton's equations), and average them, obtaining
$$
\frac{d}{dt}\langle x\rangle=\frac{\langle p\rangle}{m},
\frac{d}{dt}\langle x\rangle=-\langle \frac{dV(x)}{dx}\rangle.
$$
If we were dealing with a Harmonic oscillator, $V(x)=\frac{m\omega^2x^2}{2}$ this would be the end of the story, even without assumption about the narrow wave packet:
$$
\frac{d}{dt}\langle x\rangle=\frac{\langle p\rangle}{m},
\frac{d}{dt}\langle x\rangle=-m\omega^2\langle x\rangle.
$$
For a general case we need to do some more work. Let us assume for simplicity, that at the time moment of interest the distribution of $x$ is a Gaussian wave packet:
$$
w(x)=|\psi(x)|^2=\frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{(x-\mu)^2}{2\sigma^2}},
$$
where, of course, $\mu=\langle x\rangle$.
We can now try to evaluate the average of the force:
$$
\langle F(x)\rangle = -\langle \frac{dV(x)}{dx}\rangle=\int dx w(x) F(x)=\\
\langle \sum_{n=0}^\infty F^{(n)}(\mu)(x-\mu)^n\rangle=
\sum_{n=0}^\infty F^{(n)}(\mu)\langle (x-\mu)^n\rangle=
F(\mu) + \sum_{n=1}^\infty F^{(n)}(\mu)\langle (x-\mu)^n\rangle
$$
We note that the averages $\langle (x-\mu)^n\rangle$ for odd $n$ vanish due to the symmetry of the integral, whereas the even ones are easily calculated:
$$
\langle (x-\mu)^{2n}\rangle = \sqrt{\pi}(\pi\sigma^2)^n(2n+1)!!,
$$ (where $(2n +1)!!=1\cdot 3\cdot 5\cdot...\cdot(2n+1)$)
We thus have
$$
\langle F(x)\rangle=
F(\mu) + \sum_{n=1}^\infty F^{(2n)}(\mu)\sqrt{\pi}(\pi\sigma^2)^n(2n+1)!!
$$
Provided that the potential/force changes slowly, compared to the size of the wave packet, i.e.,
$$
\left|\frac{d^{2n}}{d\mu^{2n}}F(\mu)\right|\sigma^2\ll \frac{1}{(2n+1)!!},$$
we can neglect all the terms except the first, and we have approximately
$$
\langle F(x)\rangle=
F(\langle x\rangle).$$
This condition is quite easy to satisfy for any finite power potential, like $V(x)=kx^{2n}$, whereas for some potentials the condition breaks for high $n$, due to the rapid growth of the factorial. One can still argue that for a sufficiently narrow packet the dynamics remains nearly classical for long periods of time.
