It was suggested that this might be a better place for this question than Mathematics Stack Exchange where I asked it before.
Suppose one has a black-box function which can be evaluated anywhere (cheaply) on a specified interval $[a,b]$ and has no noise (except floating point granularity, say). What would be the best way to find the discontinuities of this function? I don't know how many discontinuities there might be and there may be none.
I can think of some straightforward methods (uniform sampling, refine where there are large differences between samples, ...), but perhaps there is a better way?
The function is "reasonable" in that one could assume that it has at most finitely many discontinuities, the same for higher derivatives, I don't mind if small pathological discontinuities are missed ... (the application is automated plotting of 1d functions).
Thanks to all who replied, particularly Pedro; the method described in Pachón, Platte and Trefethen seems to be the best approach to me, so I go now to implement it