Finding linear regimes in data

I am trying to find all linear regimes in a certain plot, given certain data points.
If you are given two vectors of data $y$ and $x$ and $y$ is some possible non-differentiable function of $x$, how would you go about finding all the linear regions in $y$? Of course, after the data has been denoised and is sufficiently clean.

The one approach I am thinking about is taking the second derivative using finite difference methods. So essentially, evaluate as frequently as possible $frac{d^2y}{dx^2} approx frac{y(n+1)-2y(n)+y(n-1)}{h^2}$ and see how close it is to $0$ upto a certain tolerance, taking into account the scale of the data. Yes, this method does not seem great, as the number of samples I have to take is rather large, and the larger the data vectors get, the longer it will take.

Could anyone suggest any papers, or alternative algorithms to determine how to find the linear regimes in a dataset?

I appreciate any advice you may have.