Which interpolation method for complicated, smooth curves? Which interpolation method should I use for complicated "smooth" curves such as 
$\frac{sin(x)}{x}$ for $x>0$.

There are many different interpolation methods and depending on the data you want to interpolate some might work better than others. If you know that the dataset comes from a trigonometric function, say $\frac{sin(x)}{x}$, then using trigonometric interpolation might be right for you.
If you don't know if the dataset comes from a specific class of functions then using splines for interpolation could work well. Here is one such method for curve interpolation.
Lets say that you have a dataset $C = (c_0, \dots, c_m)$ with $c_i \in \mathbb{R}^d$. We wish to find a B-spline
$$
f(t) = \sum\limits_{i = 0}^m p_i N_i^n(t)
$$
that interpolates $C$. This means that for a parameter vector $T = (t_0,\dots,t_m)$ we have $f(t_i) = c_i$ for all $0\leq i \leq m$. For this we need to decide on the degree $n$ of the B-spline, create the basis functions $N_i^n$ from a knot vector $U = (u_0,\dots,u_k)$ and find the control points $P = (p_0,\dots,p_m)$, $p_i \in \mathbb{R}^d$. A common choice for degree is $n = 3$ which we will use here.


*

*To create the parameter vector $T$ we will use the chord length method defined by
\begin{align*}
t_0 &= 0\\
t_i &= t_{i - 1} + \frac{|c_{i} - c_{i - 1}|}{L}\\
t_m &= 1
\end{align*}
where $L = \sum_i |c_i - c_{i - 1}|$ is the total length of the chord. The uniform or centripetal method are other parameter selection methods and a description and good discussion on their pros and cons are written by (T. Foley and G. Nielson, Knot selection for parametric spline interpolation).

*Use knot averaging suggested by de Boor (de Boor, A practical guide to splines p. 219) defined by
\begin{align*}
u_0 &= \dots = u_3 = 0\\
u_{3 + j} &= \frac{1}{3} \sum\limits_{i = j}^{3 + j - 1} t_j\\
u_{k - 3} &= \dots = u_k = 1.
\end{align*}
Here, $k$ has to satisfy $k = n + m + 1$ which means that for $n = 3$ we have $k = m + 4$. This will create a clamped B-spline.

*Create the basis functions $N_i^3$ from the knot vector $U$ and use them to calculate the matrix
\begin{align*}
A_T = \begin{pmatrix}
{N_0^3}(t_0) & {N_1^3}(t_0) & \dots & {N_m^3}(t_0) \\
{N_0^3}(t_1) & {N_1^3}(t_1) & \dots & {N_m^3}(t_1) \\
\vdots & \vdots & \ddots & \vdots \\
{N_0^3}(t_m) & {N_1^3}(t_m) & \dots & {N_m^3}(t_m)
\end{pmatrix}.
\end{align*}
Solve the linear equation system $A_T P^\top = C^\top$ (where $P^\top$ is the transpose of $P$) to get the control points $P$.

*Create the B-spline $f$ from $P$ and the basis functions $N_i^3$. This is a B-spline curve that will interpolate the dataset $C = (c_0,\dots,c_m$).
For this method to work well it is important to use a good parameter selection method. The chord length method is good for datasets that behave fairly well but others might work better for you. The paper by T. Foley and G. Nielson has a deeper discussion on the behaviour on different parameter selection methods. 
