I would like to interpolate a set of points in the real plane $(x_i,y_i), \ 1\leq i \leq n$ with specified end derivatives up to the second order. That is finding $f \in \mathcal{C}^2(\mathbb{R},\mathbb{R})$ such that $$ \forall 1 \leq i \leq n, \ f(x_i)=y_i $$ $$ f'(x_1^+) = a, \ f'(x_n^-) = b$$ $$ f''(x_1^+) = c, \ f''(x_n^-) = d$$
Using cubic splines I can either fix first or second order end derivatives but not both. Is there another family of interpolation functions that would solve my problem ?
Splines of all degrees exist and are fairly easy to compute (just by solving linear equations). The simplest interpolant would just be a polynomial -- what CAD people call a Bezier curve.
You have $n+4$ constraints, (the $n$ points plus four end conditions), so you will need a polynomial thta has $n+4$ coefficients, which means it will have degree $n+3$. Call it $f$, and call its coefficients $p_1, p_2, \ldots , p_{n+4}$. The equations you wrote give you $n+4$ linear equations that you can solve for $p_1, p_2, \ldots , p_{n+4}$.
There are some numerical difficulties, so, if you write code to do this, you need to be a bit careful.