\(\(\(
\)\(<< 
      Statistics`DescriptiveStatistics`

 (*\ \ \ \ Hermite\ spline\ interpolation\ of\ f_t\ \((x)\)\  = \ \((x - 
                t)\) _ + \(\(^\)\((k - 
                1)\)\)\ \ \ \ \ \ \ *) 
 (*\ \ \ \ on\ [0, 
            1]\ with\ 2  k - 
        4\ interior\ knots\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ *) 
 (*\ \ \ \ generated\ as\ order\ statistics\ of\ 2  
          k - 4\ i . i . d . \ uniform \((0, 
            1)\)\ \ \ \ *) 
 (*\ \ \ \ random\ variables\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ *) 
 (*\ \ \ \ Restriction : \ \ k\  ≥ \ 
          3\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ *) 


    Print["k (2k-1) degree spline interpolation"]

    k = 4

    Print["d = number of digits of accuracy"]

    d = 30



    f[x_, j_] := x^j

    g[x_, y_, 
        k_] := \((\((x - y)\)*UnitStep[x - y])\)^\((2*k - 
            1)\)

    fp[x_, j_] := j*x^\((j - 1)\)

    gp[x_, y_, 
        k_] := \((2*k - 
            1)\)*\((\((x - y)\)*UnitStep[x - y])\)^\((2*k - 
              2)\)

    h[x_, t_, 
        k_] := \((\((x - t)\)*UnitStep[x - t])\)^\((k - 1)\)/\(\((k - 
              1)\)!\)

    hp[x_, t_, 
        k_] := \((\((x - t)\)*UnitStep[x - t])\)^\((k - 2)\)/\(\((k - 
              2)\)!\)

    Basis[x_, vec_List, k_] := 
      Join[{1}, Table[f[x, j], {j, 1, 2*k - 1}], 
        Table[g[x, vec[\([j]\)], k], {j, 2, 2*k - 3}]]

    \(HermiteInterp[x_, coeff_List, vec_List, k_] := 
        coeff . Basis[x, vec, k];\)

    mxlist[coeff_List, vec_List, t_, \ k_] := 
      Table[
        Abs[HermiteInterp[i/1000, coeff, vec, k] - h[i/1000, \ t, \ k]], {i, 
          1, 1000}]


    Print["Evaluation Point t"]

    tt = Random[Real, {0, 1}, d]

    Print["interior knots"]

    VectauUnif = 
      Sort[Table[
          Random[Real, {0, 1}, d], {j, 1, 2*k - 4}]]

    Print["all knots"]

    Vectau = Append[Prepend[VectauUnif, 0], 1]

    Print["Matrix MM"]

    MM = Join[
Table[
          Join[{1}, 
Table[
              f[Vectau[\([j]\)], i], {i, 1, 2*k - 1}], 
            Table[g[Vectau[\([j]\)], Vectau[\([i]\)], k], {i, 2, 
                2*k - 3}]], {j, 1, 2*k - 2}], 
Table[
          Join[{0, 1}, 
Table[
              fp[Vectau[\([j]\)], i], {i, 2, 2*k - 1}], 
            Table[gp[Vectau[\([j]\)], Vectau[\([i]\)], k], {i, 2, 
                2*k - 3}]], {j, 1, 2*k - 2}]]

    Print["h at knot points"]

    VecH = 
      Join[Table[h[Vectau[\([j]\)], tt, k], {j, 1, 2*k - 2}], 
        Table[hp[Vectau[\([j]\)], tt, k], {j, 1, 
            2*k - 2}]]


    Print["hermite coefficients"]

    HermiteCoeff = 
      LinearSolve[MM, 
        VecH]

 (*
Print["interpolation error at i/1000, i=1,..., 1000"]
*) 

    \(mtl = mxlist[HermiteCoeff, Vectau, tt, \ k];\)

    Print["Maximum interpolation error"]

    mx = Max[mtl]

    Print["plot of function being interpolated, f_t (x) = (x-t)_+^(k-1)"]

    P1\  = \ 
      Plot[h[x, tt, k], {x, 0, 1}, PlotRange -> All]

    Print["plot of interpolant"]

    P2 = Plot[HermiteInterp[x, HermiteCoeff, Vectau, k], {x, 0, 1}, 
        PlotRange -> All]

    Print["plot of function and interpolant"]

    Show[P1, P2]

    Print["plot of interpolation error"]

    P3 = Plot[
        HermiteInterp[x, HermiteCoeff, Vectau, k] - \ h[x, tt, k], {x, 0, 
          1}, \ PlotRange -> 
          All]



    \)\)\)