Approximation with stability constant 1

I am looking for an approximation method for a typical 1-dimensional approximation problem, that has stability constant 1.

The setting

We have

• an interval $[a,b]$
• an unknown function $f$ on $[a,b]$ that is sufficiently smooth
• a partitioning of the interval $a=x_0\le x_1 \le \dots \le x_n=b$.
• known values of $f$ at knots: $y_i=f(x_i)\ (i=0\dots n)$
• Use the notation $\overline{x}=(x_0,x_1,\dots,x_n)$ and $\overline{y}=(y_0,y_1,\dots,y_n)$

We want

• An approximating function $F_{\overline{x},\overline{y}}(x)$ on interval $[a,b]$
• $F_{\overline{x},\overline{y}}(x)$ as a function of $x$ can be given in closed form, this closed form is as simple as possible (eg a spline) and can be calculated efficiently from inputs $\overline{x},\overline{y}$
• The approximation method has stability constant 1: $$\|\overline{y}-\overline{y}^\star\|_\infty\le\epsilon \ \Longrightarrow \|F_{\overline{x},\overline{y}}-F_{\overline{x},\overline{y}^\star}\|_\infty\le\epsilon$$ where the latter norm is taken on interval $[a,b]$
• $F_{\overline{x},\overline{y}}(x)$ approximates $f(x)$ to as high order of precision as possible

What I know

• Linear interpolation does this to second-order precision, so what I am looking for is a method that has more than second order precision.
• Natural cubic spline approximates to fourth-order precision but the stability condition is absolutely not satisfied even for uniform partitions.