Off[General::spell1]; 
Clear[Differences,XY0,n];
Differences[XY0_] := 
Block[{k},
  XY = XY0;
  n = Length[XY]-1;
  X = Transpose[XY][[1]];
  Y = Transpose[XY][[2]];
  h = Table[0,{n}];
  d = Table[0,{n}];
  z = Table[0,{n+1}];
  a = Table[0,{n-1}];
  b = Table[0,{n-1}];
  c = Table[0,{n-1}];
  v = Table[0,{n-1}];
  s = Table[0,{n},{4}];
  h[[1]] = X[[2]] - X[[1]];
  d[[1]] = (Y[[2]]-Y[[1]])/h[[1]];
  For[k=2,k<=n,k++,
    h[[k]] = X[[k+1]] - X[[k]];
    d[[k]] = (Y[[k+1]]-Y[[k]])/h[[k]];
    a[[k-1]] = h[[k]];
    b[[k-1]] = 2 (h[[k-1]] + h[[k]]);
    c[[k-1]] = h[[k]];
    v[[k-1]] = 6(d[[k]]-d[[k-1]])];];

Clear[TriDiagonal,a0,b0,c0,v0];
TriDiagonal[a0_,b0_,c0_,v0_] :=
Module[{k,t},
  a = a0; b=b0; c=c0; v=v0;
  For[k=2,k<=n-1,k++,
    t = a[[k-1]]/b[[k-1]];
    b[[k]] = b[[k]] - t c[[k-1]];
    v[[k]] = v[[k]] - t v[[k-1]];];
  z[[n]] = v[[n-1]]/b[[n-1]];
  For[k=n-2,1<=k,k--,
    z[[k+1]] = (v[[k]] - 
       c[[k]] z[[k+2]])/b[[k]];];];

Clear[Coeffs];
 Coeffs:=
 Block[{},
  z[[1]]  = 0;
  z[[n+1]] = 0, ];

Clear[ComputeCoeff,d,h,z,y];
ComputeCoeff[d_,h_,z_,y_] :=
Block[{k},
  For[k=1,k<=n,k++,
    s[[k,1]] = Y[[k]];
    s[[k,2]] = d[[k]]-h[[k]](2 z[[k]]+
      z[[k+1]])/6;
    s[[k,3]] = z[[k]]/2;
    s[[k,4]] = (z[[k+1]] - z[[k]])/
       (6 h[[k]]) ]; ];

Clear[Sm, SP]
S1[x_] := If[x< X[[2]],
  ((s[[1,4]](x - X[[1]])+s[[1,3]])*
    (x-X[[1]])+s[[1,2]])(x - X[[1]])+
     s[[1,1]],0] ;
Sm[m_,x_] :=
  If[X[[m]]<=x && x<X[[m+1]],
   ((s[[m,4]](x - X[[m]])+s[[m,3]])*
   (x-X[[m]])+s[[m,2]])(x - X[[m]])+
    s[[m,1]],  0] ; 
Sn[n_,x_] := If[X[[n]]<=x,
    ((s[[n,4]](x - X[[n]])+s[[n,3]])*
    (x-X[[n]])+s[[n,2]])(x - X[[n]])+
     s[[n,1]],0] ;
SP[x_] := S1[x] +
    Sum[Sm[i,x], {i,2,n-1}] + Sn[n,x];

Clear[MakeSpline];
MakeSpline :=
Block[{},
  Differences[XY];
  TriDiagonal[a,b,c,v];
  ComputeCoeff[d,h,z,Y] ];
On[General::spell1];

XY = {{0.,1.},{1.,0.},{2.,0.}, {3.,1.},{4.,2.}, {5.,2.},{6.,1.}}

MakeSpline;
dots = ListPlot[XY,PlotStyle->{PointSize[0.02]},DisplayFunction->Identity];
gr = Plot[SP[x],{x,0,6},DisplayFunction->Identity];
txt = StringForm["Natural Spline 
 y = S[x]"];
Show[gr,dots,DisplayFunction->$DisplayFunction,
PlotLabel->txt,AxesLabel->{"x","y    "}];

Clear[parts]
parts=Plot[{Sm[1,x], Sm[6,x], Sm[4,x]}, {x,0,6}, DisplayFunction->Identity];
Show[parts,dots,DisplayFunction->$DisplayFunction]

Clear[P,x];
P[m_,x_] :=
((s[[m,4]](x - X[[m]])+s[[m,3]])*
(x-X[[m]])+s[[m,2]])(x - X[[m]])+
 s[[m,1]] ;

Print[P[1,x]]; Print[P[2,x]]; 
Print[P[3,x]]; Print[P[4,x]]; 
Print[P[5,x]]; Print[P[6,x]];

piece=Plot[{P[1,x], P[2,x],P[3,x],
 P[4,x],P[5,x],P[6,x]}, 
{x, 0,6}, PlotRange->{{0,6},{-0.2,2.2}},
DisplayFunction->Identity];
Show[piece,dots,DisplayFunction->$DisplayFunction]

XY = {{0.,0.8},{0.5,0.8},{1.,1.2}, 
{8.,3.},{11.,6.},{18.,6.}, {20.5,3.5},
{24.,3.}, {25.,0}}

MakeSpline;
dots = ListPlot[XY,PlotStyle->{PointSize[0.02]}, DisplayFunction->Identity];
gr = Plot[SP[x],{x,0,26},
  PlotRange->{{0,26},{0,8.1}},  
  AspectRatio -> 0.25,
DisplayFunction->Identity];
txt = StringForm["Natural Spline 
   y = S[x]"];
Show[gr,dots,DisplayFunction->$DisplayFunction,
PlotLabel->txt,AxesLabel->{"x","y"}];

Clear[XY,dots];
f[x_] = 1 /(1+ 25x^2);
XY= Table[{0.2 i, f[0.2 i]},{i, -5, 5}]

MakeSpline;
dots = ListPlot[XY,PlotStyle->{PointSize[0.02]},DisplayFunction->Identity];
F = Plot[f[x], {x, -1,1},DisplayFunction->Identity]; 
gr = Plot[SP[x],{x,-1,1},DisplayFunction->Identity];
txt = StringForm["Natural Spline  
y = S[x]"];
Show[ gr,dots,DisplayFunction->$DisplayFunction,
PlotLabel->txt,AxesLabel->{"x","y"}];

Show[F, gr,dots,DisplayFunction->$DisplayFunction,
PlotLabel->txt,AxesLabel->{"x","y"}];