Where did the idea of hermite interpolation came from? I am given the Hermite interpolation formula directly in my text book without ANY explanations about how it was first made (obviously it was somehow constructed for the first time with some sort of intuition ) .
the formula for n+1 data from $x_0$ till $x_n$ with $f(x_0)$ till $f(x_n)$
and with $ f^{\prime}(x_0)$ till $f^{\prime}(x_n)$
$$H_{2n+1}(x) = \sum_{j=0}^n f(x_j)H_{n,j}(x) + \sum_{j=0}^n f^{\prime}(x_j)\hat H_{n,j}(x)$$
where 
$$H_{n,j} = [1 − 2(x − x_j)L^{\prime}_{n,j}(x_j)]L_{n,j}^2(x) $$
$$ \hat H_{n,j}(x) = (x-x_j) L_{n,j}^2(x) $$
I DO understand the proof and why the polynomial agrees with data and their derivatives.  
i DO understand the intuition behind Lagrange polynomials.
so I am looking for the intuition behind the formula (how it was made) specially the construction of $H$ and $\hat H$. so instead of memorizing it i can learn it!
 A: Think Lagrange interpolation and Taylor's formula.  In fact, the best way is to look at linear interpolation as the most basic Hermite (and second most basic Lagrange) and make the connection.  Suppose you want to approximate a differentiable function $f:\mathbb R\to\mathbb R$ and you give yourself two nodes: $x_0$ and $x_1$, with $x_1\neq x_0$.  
The main idea of Hermite is to take $x_1\to x_0$, but for this to succeed you cannot use the standard Lagrange basis $L_0(x):=(x-x_1)/(x_0-x_1)$ and $L_1(x):=(x-x_0)/(x_1-x_0)$ because the denominators will go to zero while the numerators stay finite.
So the trick is to replace Lagrange basis with another one (there's nothing really special about Lagrange basis except that it makes collocation easy). For example take $M_0(x):\equiv 1$ and $M_1(x):=(x-x_0)$ (which happens to be the "naive" Vandermonde basis when $x_0=0$ but consider that an accident).  Look now for two coefficients, say $c_0,c_1$ such that the linear function $c_0M_0+c_1M_1$ interpolates $f$ at $x_0,x_1$, you solve a linear system and you'll get 
$$ 
c_0=f(x_0)
\text{ and }
c_1=\frac{f(x_1)-f(x_0)}{x_1-x_0}.
$$
(Can you smell Newton's divided differences? Never mind, focus on Hermite for now.)
So in fact, posing $h=x_1-x_0$ we see that the only quantity that depends on $h$ is $c_1=\hat c_1(h)$, while $c_0,M_0,M_1$ are $h$-independent, so we have an $h$-dependent Lagrange interpolant $L^{(h)}(x)=c_0M_0(x)+\hat c_1(h)M_1(x)$.  
This is useful.
Drums rolling, we let $x_1\to x_0$ (or $h\to0$) and track the quantities that depend on $h$, i.e., $\hat c_1(h)$ which will simply converge to $f'(x_0)$, hence
the linear interpolant becomes
$$
L^{(0)}(x)=f(x_0)M_0(x)+f'(x_0)\underbrace{(x-x_0)M_0(x)}_{M_1(x)}.
$$
But this is Taylor's expansion of order 1 at $x_0$ and it is exactly what Monsieur Hermite does!
A: Both kinds of interpolation formulas rely on the superposition principle (the sum of the effect of individuals causes is the effect of the sum of the causes), and achieve a decomposition such that every point brings its own contribution. Actually, you form a basis of polynomials and linear combinations thereof.
In the case of Lagrange, consider the special case $f(x_i)=\delta_{ij}$: all ordinates but the $j^{th}$ are zero, and the latter is one. This is easily achieved by forming the product of $(x-x_i)$ and normalizing to one at $x_j$. From these $n$ basis polynomial, you can construct the interpolant for any ordinates.
The generalization to Hermite follows the same idea. You will form two families of polynomials: the first family carries the ordinates ($f(x_i)=\delta_{ij}, f'(x_i)=0$), and the second one carries the derivatives ($f(x_i)=0,f'(x_i)=\delta_{ij}$).
The rest is technical trickery, based on the idea that by squaring a Lagrange polynomial, the simple roots become double roots and the derivative vanishes at the roots, preparing candidates for the first and second family.

More precisely, $L^2_j$ achieves $f(x_i)=\delta_{ij}$, and $f'(x_i)=0$, except at $x_j$.
Let us introduce the polynomial $Z_j=(x-x_j)L_j$, such that $Z_j(x_i)=0$ and $Z_j'(x_j)=1$.
To obtain the first family, we cancel the derivative of $L_j^2$ at $x_j$ by introducing a corrective term $-Z_j(L_j^2)'$, that derives as $-Z_j'(L_j^2)'-Z_j(L_j^2)''$, i.e. $-(L_j^2)'$ at $x_j$ and $0$ elsewhere:
$$H_j=L_j^2-Z_j(L_j^2)'=(1-2(x-x_j)L_j')L_j^2.$$
To obtain the second family, we cancel $f(x_j)$ using the product $Z_jL_j$, and we have $(Z_jL_j)'=Z_j'L_j+Z_jL_j'=\delta_{ij}$ as desired. Hence:
$$\hat H_j=Z_jL_j=(x-x_j)L_j^2.$$
A: Let's call $H_{2n+1}(x) \;,\: P(x)$ .
Assume we are given a set of data $S = \{x_0 ,x_1 , \ldots , x_n \}.$
We want $P(x)$ to pass through $f(x_j)$'s and have the same derivative as $f(x)$ for all data in $S$ .
Note that, since we have $2(n+1)$ conditions for $P(x), P(x)$ should be at least of degree $2n+1$ since a $2n+1-$ degree polynomial has $2n+2$ coefficients which can be modified to meet our conditions. 
So, just like Lagrange polynomials intuition, we construct a formula like this
$$P(x) = \sum_{j=0}^n f(x_j)A_j(x) + \sum_{j=0}^n f^{\prime}(x_j) B_j(x)$$
Note1: $A_j$ corresponds to the term in the sum that it's coefficient is $f(x_i)$
and $B_j$ corresponds to the term in the sum that it's coefficient is $f^{\prime}(x_i)$.
Note2: in all formulas $x_i \in S$ and $x_j \in S$ .
Now, since we want $\forall x_i P(x_i) = f(x_i)$ then
$$1.\quad A_j(x_i) =
\begin{cases}
1,  &  {\text{ if }x_i = x_j}  \\
0, & {\text{ else } (\text{ if }x_i \neq x_j)}\\
\end{cases}$$
So, when we give $P(x)$ data $x_i,$ the corresponding term containing $f(x_i)$ in $P(x)$ appears and all the other terms of the first sum become zero.
When calculating $P(x_i),$ we don't want any of the drivatives to appear in the result so
$$2.\quad \forall x_i B_j(x_i) = 0 $$
Since $P^{\prime}(x) = \sum_{j=0}^n f(x_j)A^{\prime}_j(x) + \sum_{j=0}^n f^{\prime}(x_j) B^{\prime}_j(x),$
just as above, we determine $A^{\prime}_j$ and $B^{\prime}_j$
$$3.\quad \forall x_i A^{\prime}_j(x_i) = 0 $$
$$4.\quad B^{\prime}_j(x_i) =
\begin{cases}
1,  &  {\text{ if } x_i = x_j}  \\
0, & {\text{ else }( \text{ if }x_i \neq x_j)}\\
\end{cases}$$
 
Now, note that the behavior of $A_j$ and $B_j$ is close to that of a Lagrange basis polynomial, that is : assume $A_j = L_j$  . Then $L_j$ agrees with condition 1. but it obviously does not meet condition 3 . (also remember that $P(x)$ is of degree $2n+1$).
So, assuming each $A_j$ has the same degree, and assuming that we want to use the property of a Lagrange basis polynomial in construction of the $A_j,$ then
$$A_j = (a_j\cdot x+b_j)\cdot L^2_j$$
Note, by using $L^2$ instead of $L,$ the degree of $A$ gets from $n$ to $2n$ and by multiplying $A$ with line $(a_i\cdot x+b_i), A$ has the power $2n+1.$
We have divided $A$ into two factors:  one of degree $1$ and one of degree $2n$ .
WHY?
Since we have only conditions $1$ and $3$ for the construction of $A$ and so dividing $A$ into the predetermined(meaning it does not have unknowns like $a_j$ and $b_j$) Lagrange basis polynomial and a line $(a_j*x+b_j),$ with two unknowns should give us a unique answer for $A_j$.
This is also the case for $B_j$.
Now, assuming $A_j$ is of the form $A_j = (a_j\cdot x+b_j)\cdot L^2_j,$ and solving the equations for the conditions $1$ and $3,$ we get :
for condition$ 1$ we know
$$A_j = (a_j\cdot x+b_j)\cdot L^2_j$$
if $x_i \neq x_j$ then $A(x_i) = 0$
Since then $L_j(x_i) = 0,$ which makes $A(x_i) = 0$ regardless of $a_j$ or $b_j,$
which gives us no information about $a_j$ or $b_j$.
if $x_i = x_j$ then $A(x_j) = 1$ so
$$ a_j (x_j L_j^2(x_j)) + b_j L_j^2(x_j) = 1 $$
since $L_j(x_j) = 1$ then
$$a_j x_j + b_j = 1\\
so \; b_j = 1-a_j x_j$$
now for condition 3. we have
$$A^{\prime}_j(x) = a_j L^2_j(x) + (a_j x + b_j) 2L_j(x)L^{\prime}_j(x)$$
if $x_i \neq x_j$ again since $L_j(x_i) = 0$ we do not get any information about $a_j$ or $b_j$.
if $x_i = x_j$ then $L_j(x_j) = 1$ so
$$a_j + (a_j x_j + b_j) 2 L^{\prime}(x_j) = 0\\
   \text{by substituting} \; b_j \; \text{from above we have} \\
    a_j = -2L^{\prime}_j(x_j), $$
thus
$$ b_j = 1+2x_j L^{\prime}_j(x_j) $$
by substituting $a_j$ and $b_j$ in $A_j$ we get:
$$A_j(x) = (-2L^{\prime}_j(x_j) x + 1+2x_j L^{\prime}_j(x_j) ) L^2_j(x)$$
so,
$$A_j(x) = [1 − 2(x − x_j)L^{\prime}_{n,j}(x_j)]L_{n,j}^2(x)$$
So, we have found $A_j(x)$ which is technically the same as $H_{j}$ or $H_{n,j}.$
In the same way, using conditions $2$ and $4$ we can find  $B_j(x):$
$$ B_j(x) = (x-x_j) L_{n,j}^2(x)$$
So, we have found $\hat H_{n,j}(x)$
