This question regards the issue of interpolation of one dimension real functions.

If one has a finite set of function values and its corresponding derivatives, one could find unique continuous piecewise cubic function that interpolates those points and has the desired derivatives. This function is at least once-differentiable and in the general case is not twice-differentiable as the second derivative is discontinuous.

On the other hand a cspline-interpolation (with some boundary conditions) of those same points gives a three twice-differentiable function (but does not respect the values of the derivatives at the points).

In other words, the cspline is the only cubic spline that is twice differentiable, which makes me think that if I want a twice-differentiable interpolation *with derivatives* I will need a higher order interpolation. For example quartic or quintic.

The question is, **does interpolation with derivatives that is also twice-differentialy exists already? In practice, should it be quartic or quintic? Is there any reference or numerical implementation for it?** 

(I experimented with quartic interpolation --for values only-- in the past, but as it is well-know, even-number order interpolations produce unstable oscillations making it unsuitable for numerical approximation, which discourages me from trying quartic order interpolation)

Here it is an illustration of function interpolation of values (on the right) and derivative values (on the left)

![cubicvscspline](https://i.stack.imgur.com/02IlQ.png)

(in this example the value of the derivative is choosen to be 1 in the sample data for illustration only)