Cubic spline interpolation for the set $$\mathcal{D}=\left\{\left(0,21\right),\left(1,24\right),\left(2,24\right),\left(3,18\right),\left(4,16\right)\right\}\,.$$ All splines of the interpolation function $$f(x) = \begin{cases}-0.30357 x^3 + 3.3036x + 21& \text{if } x \in [0,1]\\-1.4821x^3 + 3.5357x^2 -0.23214 x + 22.179& \text{if } x \in (1,2]\\3.2321x^3 -24.750x^2 + 56.339 x -15.536& \text{if } x \in (2,3]\\-1.4464x^3 + 17.357 x^2 -69.982 x + 110.79& \text{if } x \in (3,4]\end{cases}
$$ are plotted and colored in black within the interval where they compose the interpolation function $f(x)$. The applied boundary condition is called “natural”.

LaTeX source

\documentclass[11pt]{article}
\usepackage[paperheight=8.5in,paperwidth=13.0in,margin=0in,headheight=0.0in,footskip=0in,includehead,includefoot]{geometry}

% graphics
\usepackage{tikz}
\usepackage{pgfplots}
\pgfplotsset{width=10cm,compat=1.12}
\usepgfplotslibrary{fillbetween}

\begin{document}
\begin{minipage}{.35\textwidth}	\begin{tikzpicture}
		\pgfplotsset{
			scale only axis,
			xmin=-.5,xmax=4.5,
			ymin=15,ymax=25.5,
		}

		\begin{axis}[
			xlabel=$x$,
			ylabel=$y$,
			samples=100,
			xtick={0,1,2,3,4}
			]\addplot [only marks] table {
				0 21
				1 24
				2 24
				3 18
				4 16
				};
			\addplot[color=gray][]{+-0.30357142857142855*x^3+3.3035714285714284*x^1+21*x^0};
			\addplot[][domain=0:1]{+-0.30357142857142855*x^3+3.0785714285714284e-61*x^2+3.3035714285714284*x^1+21*x^0};
		\end{axis}
	\end{tikzpicture}
\end{minipage}
\begin{minipage}{.35\textwidth}
	\begin{tikzpicture}
		\pgfplotsset{
			scale only axis,
			xmin=-.5,xmax=4.5,
			ymin=15,ymax=25.5,
		}

		\begin{axis}[
			xlabel=$x$,
			ylabel=$y$,
			samples=100,
			xtick={0,1,2,3,4}
			]\addplot [only marks] table {
				0 21
				1 24
				2 24
				3 18
				4 16
				};
			\addplot[color=gray][domain=-2:3]{+-1.4821428571428572*x^3+3.5357142857142856*x^2+-0.23214285714285715*x^1+22.178571428571427*x^0};
			\addplot[][domain=1:2]{+-1.4821428571428572*x^3+3.5357142857142856*x^2+-0.23214285714285715*x^1+22.178571428571427*x^0};
		\end{axis}
	\end{tikzpicture}
\end{minipage}\\\\

\begin{minipage}{.35\textwidth}
	\begin{tikzpicture}
		\pgfplotsset{
			scale only axis,
			xmin=-.5,xmax=4.5,
			ymin=15,ymax=25.5,
		}
		\begin{axis}[
			xlabel=$x$,
			ylabel=$y$,
			samples=100,
			xtick={0,1,2,3,4}
			]\addplot [only marks] table {
				0 21
				1 24
				2 24
				3 18
				4 16
				};
			\addplot[color=gray][domain=0:4.5]{+3.232142857142857*x^3+-24.75*x^2+56.339285714285715*x^1+-15.535714285714286*x^0};
			\addplot[][domain=2:3]{+3.232142857142857*x^3+-24.75*x^2+56.339285714285715*x^1+-15.535714285714286*x^0};
		\end{axis}
	\end{tikzpicture}
\end{minipage}
\begin{minipage}{.35\textwidth}
	\begin{tikzpicture}
		\pgfplotsset{
			scale only axis,
			xmin=-.5,xmax=4.5,
			ymin=15,ymax=25.5,
		}
		\begin{axis}[
			xlabel=$x$,
			ylabel=$y$,
			samples=100,
			xtick={0,1,2,3,4}
			]\addplot [only marks] table {
				0 21
				1 24
				2 24
				3 18
				4 16
				};
			\addplot[color=gray][domain=2:5]{+-1.4464285714285714*x^3+17.357142857142858*x^2+-69.98214285714286*x^1+110.78571428571429*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{minipage}
\end{document}