Interpolation and divergence of $n$-th derivative i have a curiosity.
Let's say i have a function $f \in C^{\infty}[a,b]$ such that there's a $x_0 in [a,b]$ that makes and $f(x_0) = 0$ and $a_n = f^{(n)}(x_0)$ diverges (i.e. $|a_n|=\infty$) could this affect somehow the accuracy of the approximation if i interpolate a function?
 A: The error of interpolation is usually expressed in terms of the higher derivatives of the interpolated function, similar to Lagrange  Remainder in terminated Taylor expansion which also sort of interpolation, actually extrapolation.
Namely, if $I_n(x)$ is an $n$'th order interpolant of $f(x)$ at points $\{(x_i,f(x_i)\}_{i=0}^n$  then 
$$e_n(x)=f(x)-I_n(x)\approx\frac{f^{(n+1)}(c)}{(n+1)!}\prod_{i=0}^n(x-x_i)
$$
Therefore if the derivative $|f^{(n+1)}|$ doesn't exists or unbounded all the process of interpolation isn't valid anymore as the error term isn't defined or unbounded.
If all derivatives of order  $\leq n+1$ exists and bounded then your interpolation will work, otherwise - not. 
Note that if you use too many points even if your function all required derivatives this won't be too good idea as the accuracy depends on a polynomial of a degree $n$ which grow up fast. Therefore practically one uses piecewise interpolations one per each few points (usually 2-3) of data,  or alternatively splines, the former can be considered the non smooth version of the later.
Therefore the smoothness requirements on $f$ are not too problematic.

example for $f=3^x-1$
>> x=linspace(-1,1,100); % uniform grid interpolation points
>> y=linspace(-1,1,2000);
>> f=@(x) 3.^x-1
>>
>>g=interp1(x,f(x),y,'linear');
>> norm(f(y)-g)

ans =

   0.002859258860134

>> g=interp1(x,f(x),y,'spline');
>> norm(f(y)-g)

ans =

     7.176506074372971e-08

>> g=interp1(x,f(x),y,'pchip');
>> norm(f(y)-g)

ans =

     5.546115531849058e-06

switching to chebyshev points (linear interpolation has a problem at end points)
>> x=cos((2*(1:100)-1)*pi/200);
>> g=interp1(x,f(x),y,'linear');
>> norm(f(y)-g)

ans =

   NaN

>> norm(f(y(2:end-1))-g(2:end-1))

ans =

   0.004192591427323


>> g=interp1(x,f(x),y,'spline');
>> norm(f(y)-g)

ans =

     7.481779555065644e-08

>> g=interp1(x,f(x),y,'pchip');
>> norm(f(y)-g)

ans =

     1.639113701734916e-05

