Procedimento

> restart;
with(linalg):

Warning, new definition for norm

Warning, new definition for trace

> CubSpline4Ptsf:=proc(X::vector,f::procedure)
local h,M,fX,A,b,sol,U,S1,S2,S3,spline,i,S,dim;
dim:=vectdim(X);
h:=vector(dim-1,k->X[k+1]-X[k]);
M:=vector(dim,0):
fX:=vector(dim,[f(X[1]),f(X[2]),f(X[3]),f(X[4])]);
A:=matrix(2,2,0):
A[1,1]:=(h[1]+h[2])/3:
A[1,2]:=h[2]/6:
A[2,1]:=h[2]/6:
A[2,2]:=(h[2]+h[3])/3:
A:=evalm(A);
b:=vector(2,0):
b[1]:=(fX[3]-fX[2])/h[2]-(fX[2]-fX[1])/h[1]:
b[2]:=(fX[4]-fX[3])/h[3]-(fX[3]-fX[2])/h[2]:
b:=evalm(b);
sol:=vector(2,0):
sol:=Gauss(A,b):
M[2]:=sol[1];
M[3]:=sol[2];

for i from 1 to 3 do
S[i]:=(X[i+1]-x)^3/(6*h[i])*M[i]+(x-X[i])^3/(6*h[i])*M[i+1]+(fX[i]/h[i]-h[i]/6*M[i])*(X[i+1]-x)+(fX[i+1]/h[i]-h[i]/6*M[i+1])*(x-X[i]);
S[i]:=simplify(S[i]);
od;
spline:=piecewise(X[1]<=x and x<X[2],S[1],X[2]<=x and x<X[3],S[2],X[3]<=x and x<=X[4],S[3]);
RETURN(spline);
end:

> X:=vector(4,[-1,-1/3,1/3,1]);
dim:=vectdim(X);
h:=vector(dim-1,k->X[k+1]-X[k]);

> f:=x->1/(1+25*x^2);
CubSpline4Ptsf(X,f);

> p:=x->CubSpline4Ptsf(X,f):
plot({p(x),f(x)},x=X[1]..X[4]);
plot({p(x)},x=X[1]..X[4]);

>