Variante 1 -- Exercício 4.18

> h0:=2/3:h1:=2/3:h2:=2/3:h3:=2/3:
h:=vector(3,[2/3,2/3,2/3]);
M:=vector(4,0):
f:=x->1/(1+25*x^2);
X:=vector(4,[-1,-1/3,1/3,1]);
fX:=vector(4,[f(X[1]),f(X[2]),f(X[3]),f(X[4])]);

> A:=matrix(2,2,0):
A[1,1]:=(h0+h1)/3:
A[1,2]:=h1/6:
A[2,1]:=h1/6:
A[2,2]:=(h1+h2)/3:
A:=evalm(A);
b:=vector(2,0):
b[1]:=(fX[3]-fX[2])/h1-(fX[2]-fX[1])/h0:
b[2]:=(fX[4]-fX[3])/h2-(fX[3]-fX[2])/h1:
b:=evalm(b);

> sol:=vector(2,0):
#U:=gausselim(A):
#sol:=backsub(U,b):
sol:=Gauss(A,b);
M[2]:=sol[1];
M[3]:=sol[2];

> h:=2/3:
S1:=(x-X[1])^3/(6*h)*M[2]+(fX[1]/h)*(X[2]-x)+(fX[2]/h-h/6*M[2])*(x-X[1]);
S2:=(X[3]-x)^3/(6*h)*M[2]+(x-X[2])^3/(6*h)*M[3]+(fX[2]/h-h/6*M[2])*(X[3]-x)+(fX[3]/h-h/6*M[3])*(x-X[2]);
S3:=(X[4]-x)^3/(6*h)*M[3]+(fX[3]/h-h/6*M[3])*(X[4]-x)+fX[4]/h*(x-X[3]);

> S1:=simplify(S1);
S2:=simplify(S2);
S3:=simplify(S3);

> piecewise(X[1]<=x and x<X[2],S1,X[2]<=x and x<X[3],S2,X[3]<x and x<=X[4],S3);
p:=x->piecewise(X[1]<=x and x<X[2],S1,X[2]<=x and x<X[3],S2,X[3]<=x and x<=X[4],S3):

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

> with(student):

> fS1:=makeproc(S1,x);
DfS1:=D(fS1);
DfS1(1/3);
DfS1(-1/3);evalf(DfS1(-1/3));

> fS2:=makeproc(S2,x);
DfS2:=D(fS2);
DfS2(-1/3);
S2X2:=fS2(1/3);
dS2X2:=DfS2(1/3);evalf(DfS2(1/3));

> fS3:=makeproc(S3,x);
A:=plot(fS3(x),x=1/3..1):
B:=plot(fS2(x),x=-1/3..1/3):
with(plots):
display(A,B);
DfS3:=D(fS3);
S3X2:=fS3(1/3);
dS3X2:=DfS3(1/3);evalf(DfS3(1/3));

>