A curve pass via points at TiKz

You can use plot [smooth] coordinates (which is not a single polynom but a spline):

\documentclass[tikz]{standalone}

\begin{document}
\begin{tikzpicture}
    \draw[style=help lines] (-5,-5) grid (5,5);
    \draw (-4,0)--(4,0);
    \draw (0,-4)--(0,4);
    \foreach \y in {-4,-3,...,4} {
        \draw (0 - 0.1,\y) -- (0+0.1,\y);
        \draw (\y,0 - 0.1) -- (\y,0+0.1);
    }

    %Nodes:
    \node (a0) at (-4,-4) {};
    \draw[fill] (a0) circle [radius=1.5pt];
    \node (a1) at (-2,4) {};
    \draw[fill] (a1) circle [radius=1.5pt];
    \node (a2) at (2,-2) {};
    \draw[fill] (a2) circle [radius=1.5pt];
    \node (a3) at (4,2) {};
    \draw[fill] (a3) circle [radius=1.5pt];         

    \draw plot [smooth] coordinates {(-4,-4)  (-2,4)  (2,-2)  (4,2)}; % to (a2) to (a3);

\end{tikzpicture}
\end{document}

Solution which forces the middle points to have a horizontal tangent:

\documentclass[tikz,border=3.14]{standalone}

\begin{document}
\begin{tikzpicture}
    \draw[style=help lines] (-5,-5) grid (5,5);
    \draw (-4,0)--(4,0);
    \draw (0,-4)--(0,4);
    \foreach \y in {-4,-3,...,4} {
        \draw (0 - 0.1,\y) -- (0+0.1,\y);
        \draw (\y,0 - 0.1) -- (\y,0+0.1);
    }

    %Nodes:
    \node (a0) at (-4,-4) {};
    \draw[fill] (a0) circle [radius=1.5pt];
    \node (a1) at (-2,4) {};
    \draw[fill] (a1) circle [radius=1.5pt];
    \node (a2) at (2,-2) {};
    \draw[fill] (a2) circle [radius=1.5pt];
    \node (a3) at (4,2) {};
    \draw[fill] (a3) circle [radius=1.5pt];         

    \draw (-4,-4) to[out=90,in=180] (-2,4) to[out=0,in=180] (2,-2) to[out=0,in=-95] (4,2); % to (a2) to (a3);

\end{tikzpicture}
\end{document}

I don't know how to compute this in LaTeX easily, so I fitted a plot using Python's numpy.polyfit and used the result to plot the fit in TikZ:

\documentclass[tikz,border=3.14]{standalone}

%% polynomial coefficients found with Python (numpy.polyfit)
%% $f(x) = 0.1875 x^3 - 1/6 x^2 - 2.25 x^1 + 10/6 x^0$

\begin{document}
\begin{tikzpicture}
    \draw[style=help lines] (-5,-5) grid (5,5);
    \draw (-4,0)--(4,0);
    \draw (0,-4)--(0,4);
    \foreach \y in {-4,-3,...,4} {
        \draw (0 - 0.1,\y) -- (0+0.1,\y);
        \draw (\y,0 - 0.1) -- (\y,0+0.1);
    }

    %Nodes:
    \node (a0) at (-4,-4) {};
    \draw[fill] (a0) circle [radius=1.5pt];
    \node (a1) at (-2,4) {};
    \draw[fill] (a1) circle [radius=1.5pt];
    \node (a2) at (2,-2) {};
    \draw[fill] (a2) circle [radius=1.5pt];
    \node (a3) at (4,2) {};
    \draw[fill] (a3) circle [radius=1.5pt];         

    \draw plot[domain=-4:4,samples=100] (\x, .1875*\x*\x*\x - \x*\x/6 - 2.25*\x + 10/6);

\end{tikzpicture}
\end{document}

Just for your information. You can calculate and plot the interpolation polynomial with Python and the two libraries Matplotlib and NumPy:

import numpy as np
import matplotlib.pyplot as plt

x = (-4, -2, 2, 4)
y = (-4, 4, -2, 2)
p = np.polyfit(x,y,3)
t = np.linspace(min(x),max(x),num=100)
f = np.polyval(p,t)
plt.plot(t,f)

Matplotlib supports export to TikZ code (actually it exports to PGF) and to save the plots directly as PDF created with TikZ and LaTeX (see for example https://tex.stackexchange.com/a/426071/117050 and https://tex.stackexchange.com/a/391078/117050 for some code that might get you started).

With some calculations, I found formula of the function is -(1/72)*\x^4+3/16*(\x^3)+(1/9)*\x^2-9/4*\x+7/9 I use pgfplots to draw

\documentclass[tikz]{standalone}
\usepackage{pgfplots}
\pgfplotsset{compat=1.16}
\usepackage{fouriernc}
\begin{document}
    \begin{tikzpicture}[
    declare function={
        f(\x)=-(1/72)*\x^4+3/16*(\x^3)+(1/9)*\x^2-9/4*\x+7/9;
    }
    ]
    \begin{axis}[axis equal,
    width=12 cm, 
    grid=major, 
    axis x line=middle, axis y line=middle,
    axis line style = very thick,
    grid style={gray!30},
    ymin=-5, ymax=5, yticklabels={}, ylabel=$y$,
    xmin=-5, xmax=5, xticklabels={}, xlabel=$x$,
    samples=500,
    ]
    \addplot[blue, very thick,domain=-5:5, smooth]{f(x)};
    \node[below] at (-2, 0) {$-2$};
    \node[above ] at (-4, 0) {$-4$};
    \node[below ] at (4, 0) {$4$};
    \node[right] at (0,-4) {$-4$};
    \node[left  ] at (0,2) {$2$};
    \node[ right ] at (0,4) {$4$};
    \node[below right] at (0, 0) {$O$};
    \node[above ] at ( 2,0) {$2$};
    \node[left  ] at (0, -2) {$-2$};
    \addplot [mark=*,only marks,samples at={-4,-2,2,4}] {f(x)};
    ;
    \draw[dashed, thick] (-4,0) -- (-4,-4) -- (0,-4);
    \draw[dashed, thick] (-2,0) -- (-2,4) -- (0,4);
    \draw[dashed, thick] (2,0) -- (2,-2) -- (0,-2);
    \draw[dashed, thick] (4,0) -- (4,2) -- (0,2);
    \end{axis}
    \end{tikzpicture}
\end{document} 

Results from Maple.

With marmot's help , I reduce my code

\documentclass[tikz]{standalone} 
\usepackage{pgfplots} 
\pgfplotsset{compat=1.16} 
\usepackage{fouriernc} 
\begin{document} 
\begin{tikzpicture}[ 
declare function={ 
f(\x)=-(1/72)*\x^4+3/16*(\x^3)+(1/9)*\x^2-9/4*\x+7/9; 
} 
] 
\begin{axis}[axis equal, 
width=12 cm, 
grid=major, 
axis x line=middle, axis y line=middle, 
axis line style = very thick, 
grid style={gray!30}, 
ymin=-5, ymax=5, yticklabels={}, ylabel=$y$, 
xmin=-5, xmax=5, xticklabels={}, xlabel=$x$, 
samples=500, 
] 
\addplot[blue, very thick,domain=-5:5, smooth]{f(x)}; 

\addplot [mark=*,only marks,samples at={-4,-2,2,4}] {f(x)}; 
; 

\pgfplotsinvokeforeach{-4,-2,2,4}{\draw[dashed] ({#1},0) |- (0,{f(#1)}); } 

\foreach \X/\Y in {-4/right,-2/left,2/left,4/right} 
{\edef\temp{\noexpand\node[\Y] at (0,\X) {$\X$};} 
\temp} 


\foreach \X/\Y in {-4/above,-2/below,2/above,4/below} 
{\edef\temp{\noexpand\node[\Y] at (\X,0) {$\X$};} 
    \temp} 
% 
\end{axis} 
\end{tikzpicture} 
\end{document}

Another way

\documentclass[tikz,12pt]{standalone} 
\usepackage{pgfplots} 
\pgfplotsset{compat=1.16} 
\usepackage{fouriernc} 
\begin{document} 
\begin{tikzpicture}[ 
declare function={ 
f(\x)=-(1/72)*pow(\x,4)+3/16*pow(\x,3)+(1/9)*\x*\x-9/4*\x+7/9; 
xmin=-5;xmax=5;ymin=-5;ymax=5;} 
] 
\draw[gray!30] (xmin,ymin) grid (xmax,ymax); % grid 
\draw[->, thick] (xmin,0)--(xmax,0) node [below left]{$x$}; 
\draw[->,thick] (0,ymin)--(0,ymax) node [below left]{$y$}; 
\foreach \X in {-4,-2,2,4} {\draw[dashed] (\X,0) |- (0,{f(\X)}); } 
\node[below right] at (0, 0) {$O$}; 
\foreach \Y in {-4,-2,2,4} \fill (\Y,{f(\Y)}) circle(2pt); 
\foreach \p/\g in {-4/90,-2/-90,2/90,4/-90 }\draw(\p,0)node[shift={(\g:.3)},scale=1]{$\p$}--+(0,.05)--+(0,-.05); 
\foreach \p/\g in {-4/0,-2/180,2/45,4/0}\draw(0,\p)node[shift={(\g:.3)},scale=1]{$\p$}--+(0,.05)--+(0,-.05); 
\clip (xmin,ymin) rectangle (xmax,ymax);
\draw[smooth,samples=300,very thick, blue] plot(\x,{f(\x)});    \end{tikzpicture} 
\end{document}

We can use \draw controls - the red curve, in comparison with the blue curve \draw plot[smooth] coordinates. (if you want, you can control so that the red and blue curves are identical)

\documentclass[tikz,border=5mm]{standalone}
\begin{document}
\begin{tikzpicture}
\draw[gray!30] (-5,-5) grid (5,5);
\draw (-5,0)--(5,0) (0,-5)--(0,5);
\foreach \i in {-5,...,5}
\draw 
(0,\i)--+(1mm,0)--+(-1mm,0)
(\i,0)--+(0,1mm)--+(0,-1mm);

\draw[blue] plot[smooth] coordinates 
{(-4,-4)  (-2,4)  (2,-2)  (4,2)};

\draw[red] 
(-4,-4)..controls +(80:1) and +(180:1)..
(-2,4)..controls +(0:1) and +(180:1)..
(2,-2)..controls +(0:1) and +(-100:1)..
(4,2);

\foreach \p in {(-4,-4),(-2,4),(2,-2),(4,2)}
\fill \p circle(2pt);
\end{tikzpicture}
\end{document}