Following Christopher2222 request, I wrote the following procedures for "exact" cubic Hermite spline interpolation,

p:=proc(x0,p0,m0,x1,p1,m1,x)
local t,d;
d:=x1-x0;
t:=(x-x0)/d;
p0+(d*m0+(3*(p1-p0)-d*(2*m0+m1)+(2*(p0-p1)+d*(m0+m1))*t)*t)*t
end:

pb:=proc(x0,p0,x1,p1,m1,x)
local t,d;
d:=x1-x0;
t:=(x-x0)/d;
p0+(2*(p1-p0)-d*m1+(p0-p1+d*m1)*t)*t
end:

pe:=proc(x0,p0,m0,x1,p1,x)
local t,d;
d:=x1-x0;
t:=(x-x0)/d;
p0+(d*m0+(p1-p0-d*m0)*t)*t
end:

df:=proc(L)
local A,n,i;
n:=op([2,1,2],L);
A:=Array(1..n);
for i from 2 to n-1 do 
if (L[i+1,2]-L[i,2])*(L[i,2]-L[i-1,2])>0 then 
A[i]:=(L[i+1,2]-L[i-1,2])/(L[i+1,1]-L[i-1,1]) fi od;
A
end:

Chris:=proc(L,X)
local x,m,n,i;
m:=df(L);
n:=op([2,1,2],L);
eval(piecewise(x<L[2,1],pb(L[1,1],L[1,2],L[2,1],L[2,2],m[2],x),
seq('x<L[i+1,1],p(L[i,1],L[i,2],m[i],L[i+1,1],L[i+1,2],m[i+1],x)',i=2..n-2),
pe(L[n-1,1],L[n-1,2],m[n-1],L[n,1],L[n,2],x)),x=X)
end:

It's just a rough sketch, unfinished, not covering all the cases, and not tested for bugs.

Here is an example of its usage,

L := [3, 4, 6, 7, 2, 3, 5, 4, 6, 8, 20, 4, 5, 12, 0, 5, 5, 5, 3]:
b := Array(ListTools:-Enumerate(L)):
plots:-display(plot(Chris(b,x),x=1..19),plots:-pointplot(b));

_______________
Alec Mihailovs, PhD