Using Lagrange's interpolation formula $$\begin{equation*} P_N(x)=\sum_{i=0}^{N}\prod_{k\ne i} \frac{x-x_k}{x_i-x_k}y_i, \end{equation*}$$ we could attempt to approximate the function $f(x)$ with a first-order polynomial in $x$ in the two sub-intervals $x\in[x_0-h,x_0]$ and $x\in[x_0,x_0+h]$. A first order polynomial means simply that we have for say the interval $x\in[x_0,x_0+h]$ $$\begin{equation*} f(x)\approx P_1(x)=\frac{x-x_0}{(x_0+h)-x_0}f(x_0+h)+\frac{x-(x_0+h)}{x_0-(x_0+h)}f(x_0), \end{equation*}$$ and for the interval $x\in[x_0-h,x_0]$ $$\begin{equation*} f(x)\approx P_1(x)=\frac{x-(x_0-h)}{x_0-(x_0-h)}f(x_0)+\frac{x-x_0}{(x_0-h)-x_0}f(x_0-h). \end{equation*}$$