Example

  u={0,.1125,.2763,.3794}
  h={.1,.15,.1}

From our previous efforts (make sure you see where these came from!) , our equations are:
.2s0+.1s1=3(.1125)=.3375

10s0 + 100/3 s1 + 20/3 s2
= 3/(.15)^2 (.2763) + (3/(.15)^2 - 3/(.1)^2)(.1125) - 3/(.1)^2(0)
=18.09

1/(.15)^2 s1 + (2/(.15)^2 + 2/(.1)^2)s2 + 1/(.1)^2 s3
= 3/(.1)^2 (.3794) + (3/(.15)^2 - 3/(.1)^2)(.2763) - 3/(.15)^2 (.1125)
=52.77

.1s2+.2s3=3(.1031)=.3093

  matrix={{2/10,1/10,0,0},{10,100/3,20/3,0},
  {0,(20/3)^2,2(20/3)^2-300,100},{0,0,1/10,2/10}};
  MatrixForm[matrix]

  du={.3375,18.09,52.77,.3039}

  s=LinearSolve[matrix,du]

This gives s0, s1, s2, s3, which Mathematica sees as s[[1]], s[[2]], s[[3]], and s[[4]]

To the left of (0,0) we have a line with slope s0=1.70985 through (0,0).

  Clear[x]

  newspline0[x_]:=s[[1]]x;
  newspline0[x]

Between x=0 and x=0.1 we have

  newspline1[x_]:=genherm[x,0,s[[1]],.1125,s[[2]],h[[1]]];
  newspline1[x]

Between x=.1 and x=.25 we have

  newspline2[x_]:=gc2[x,.1125,s[[2]],.2763,s[[3]]]/.{h1->h[[1]],h2->h[[2]]};
  newspline2[x]

In our Theory section we did not go beyond gc2, so we must find gc3 (spline3) to be between x=.25 and x=.35.

  Clear[u]
  DSolve[{u''''[x]==0,u[.25]==.2763,u'[.25]==s[[3]],
  		u[.35]==.3794,u'[.35]==s[[4]]},u[x],x]/.x->x

  newspline3[x_]:=
   1.374566719314079424 - 11.38577084837545128*x + 
   
       36.88044584837545133*x^2 - 35.63851985559566794*x^3

To the right of x=.35 we have a line with slope 1.33339 through (.35, .3795)

  newspline4[x_]:=.3795+s[[4]](x-.35);
  newspline4[x]

Let's put it together!

  newspline[x_]:=If[x<=0,newspline0[x],
  				If[x<=.1,newspline1[x],
  					If[x<=.25,newspline2[x],
  						If[x<=.35,newspline3[x],newspline4[x]]
  						
  					]
  				]
  			]

  xvals={0,.1,.25,.35}
  u={0,.1125,.2763,.3794}

  plot3=ListPlot[Table[{xvals[[i]],u[[i]]},{i,1,4}],PlotStyle->{RGBColor[1,0,0],PointSize[.02]}];
  plot4=Plot[newspline[x],{x,-.1,.4}];
  Show[plot3,plot4];

The points that we have been considering actually are approximations from the error function Erf[x] which is defined by
Erf[x]=2/Sqrt[Pi] Integral from 0 to x of Exp[-t^2]dt.
Let's see how we did.

  plot5=Plot[Erf[x],{x,-.1,.4}];
  Show[plot3,plot4,plot5];

Up to Best Bet: Cubic Spline Interpolation