$f'(x)=0$ makes every other (whole) derivative zero, yet this can occur at a single point. Finite difference fails since values can change in-between.
We can tell if a function is constant over a selected interval, $x_0$ to $x_1$, via inner product with an invertible family of orthogonal basis functions - e.g. its Fourier transform will be non-zero only for $\omega=0$ (where we extend $f(x_0)$ to $-\infty$ and $f(x_1)$ to $\infty$). But can every such $(x_0, x_1)$ be found without brute-force inspection?
If necessary, assume "nice" behavior (e.g. polynomial-like); I don't seek the most general case. For clarity, I mean continuously constant over a finite interval.
One approach is normalized cross-correlation with a brickwall (pillbox/boxcar/aperiodic square pulse) at every possible width $\Delta x$ (or just those of interest); the result will peak only where $(x_0, x_0 + \Delta x)$ is constant. For one width and one interval:
$$ K(\Delta x)\left<f(x), B(x - x_0, \Delta x) \right> = K(\Delta x)\int_{-\infty}^{\infty} f(x) B(x - x_0, \Delta x)dx $$
where $K$ is a width-dependent (for $B$) product of L2 norms:
$$ K(\Delta x) = \frac{1}{||f(x)|| \cdot||B(x, \Delta x)||} $$
Repeating this for every $x_0$ and $\Delta x$ defines the "boxcar transform" (let me know if it has a name):
$$ \mathbb{BOX}_{f(x)}(\Delta x, x_0) = \int_{-\infty}^{\infty} K(\Delta x) \left( \int_{-\infty}^{\infty} \left( \int_{-\infty}^{\infty} f(x) B(x - x_0, \Delta x)dx \right) dx_0 \right) d\Delta x $$ which yields a two-dimensional representation, where for some $\Delta x$, a peak at $x_0$ indicates a constant interval $(x_0, x_0 + \Delta x)$. This can be implemented efficiently with FFT convolution.
I've not considered every edge case (e.g. $f(x)=0$), but seems this'll work in general. Alternative solutions welcome.