How to to generate Hermite interpolating polynomials? How is it possible to generate Hermite interpolating polynomial (spline) with an arbitrary degree which interpolates between $x_0$ and $x_1$ and also between the their first $((k-1)/2)$-th derivatives $m_{i,0}$ and $m_{i,1}$?
Let

*

*$k$ is odd and it denotes the degree of the polynomial

*$t \in [0,1]$ is the interpolation parameter

Then:

*

*$k=1$ is a simple linear spline
$tx_1 + (1-t)x_0$


*$k=3$ is a cubic Hermite spline
$t^3 (m_0+m_1+2 x_0-2 x_1) - t^2 (2 m_0 + m_1 + 3x_0 - 3x_1) + t m_0 + x_0$


*$k=5$ is a quintic Hermite spline
$-\frac{ 1 }{2} t ^ 5 (6 m_0 + 6 m_1 + n_0 - n_1 + 12 x_0 - 12 x_1) + \frac{ 1 }{2} t ^ 4 (16 m_0 + 14 m_1 + 3 n_0 - 2 n_1 + 30 x_0 - 30 x_1) - \frac{ 1 }{2} t ^ 3 (12 m_0 + 8 m_1 + 3 n_0 - n_1 + 20 x_0 - 20 x_1) + m_0 t + \frac{ n_0 t ^ 2 }{2}+x_0$


*and so on...
 A: There is a well-known algorithm for just polynomial interpolation, that generates polynomials in the Newton form. It goes like this. Given a function $f$ and distinct points $x_0,x_1,\dots,x_n$ we define the divided difference $f[x_i,x_{i+1},\dots,x_j]$ recursively as $\frac{f[x_{i+1},\dots,x_j]-f[x_i,\dots,x_{j-1}]}{x_j-x_i}$, with the base case $f[x_i]=f(x_i)$. Then the polynomial interpolant is 
$$f[x_0]+f[x_0,x_1](x-x_0)+\dots+f[x_0,x_1,\dots,x_n] \prod_{i=0}^{n-1} (x-x_i).$$
In hand calculations this can be conveniently evaluated using a triangular array: we make a column of $x$ values and another of $y$ values, then a column of first differences, then second differences, etc. Each column after the first two is one row shorter until we get to the last column which is of length $1$. The coefficients are then the first row of the array.
For Hermite interpolation, we duplicate an $x$ value where we want to enforce $k$ derivatives $k+1$ times in the aforementioned array. Consequently we will be asked for a divided difference that requires a division by zero. In these cases we interpret it as a limit as one $x$ value approaches another, so we substitute the given derivative value there and then continue. So for example $f[x_0,x_0]=f'(x_0)$. For another example, $f[x_0,x_0,x_1]=\frac{f[x_0,x_1]-f'(x_0)}{x_1-x_0}=\frac{f(x_1)-f(x_0)-f'(x_0)(x_1-x_0)}{(x_1-x_0)^2}$.
Note that in this framework it is not possible to solve a problem such as $p(0)=1,p''(0)=2,p(1)=3$. To handle that in this framework you would need to introduce $p'(0)$ as a free parameter and obtain a one-parameter family of cubic polynomials satisfying the desired equations. One could then choose $p'(0)$ at the end as desired (for example, to obtain a quadratic solution).
