Cubic spline interpolation with boundary condition “natural” for the set of points $$\mathcal{D}=\left\{\left(0,21\right),\left(1,24\right),\left(2,24\right),\left(3,18\right),\left(4,16\right)\right\}\,.$$The plotted function is described by $$f(x) = \begin{cases}-0.30357 \cdot x^3 + 3.3036\cdot x + 21& \text{if } x \in [0,1]\\-1.4821\cdot x^3 + 3.5357\cdot x^2 + -0.23214 \cdot x + 22.179& \text{if } x \in (1,2]\\3.2321\cdot x^3 + -24.750 \cdot x^2 + 56.339 \cdot x + -15.536& \text{if } x \in (2,3]\\-1.4464\cdot x^3 + 17.357 \cdot x^2 + -69.982 \cdot x + 110.79& \text{if } x \in (3,4].\end{cases}
$$

LaTeX source

\documentclass[11pt]{article}

% graphics
\usepackage{tikz}
\usepackage{pgfplots}
\pgfplotsset{compat=1.12}
\usepgfplotslibrary{fillbetween}


\begin{document}
	\begin{tikzpicture}

		\pgfplotsset{
			scale only axis,
		}

		\begin{axis}[
			xlabel=$x$,
			ylabel=$y$,
			samples=100,
			]\addplot [only marks] table {
0 21
1 24
2 24
3 18
4 16
};
			
			\addplot[][domain=0:1]{+-0.30357142857142855*x^3+3.0785714285714284e-61*x^2+3.3035714285714284*x^1+21*x^0};

			\addplot[][domain=1:2]{+-1.4821428571428572*x^3+3.5357142857142856*x^2+-0.23214285714285715*x^1+22.178571428571427*x^0};

			\addplot[][domain=2:3]{+3.232142857142857*x^3+-24.75*x^2+56.339285714285715*x^1+-15.535714285714286*x^0};

			\addplot[][domain=3:4]{+-1.4464285714285714*x^3+17.357142857142858*x^2+-69.98214285714286*x^1+110.78571428571429*x^0};
		\end{axis}

	\end{tikzpicture}

\end{document}