write(%io(2),'splines cubiques sur n points')
//
// profil exact
//
deff('y=f(x)','y=x.*cos(18*x)');
//
n=evstr(x_dialog('nombre de points?','15'));
timer();
n=n-2;
//
//matrice du systeme
//
v=zeros(1:n-2);
A=toeplitz([4,1,v]);
//
//points d'interpolation
//
h=1/(n+1);
x=0:h:1;
//
//definition du second membre du systeme
//
b=zeros(n,1);
c0=(f(x(2))-f(x(1)))/h;
c1=(f(x(n+2))-f(x(n+1)))/h;
b(1)=(3/h)*(f(x(3))-f(x(1)))-c0;
for i=2:n-1
	b(i)=(3/h)*(f(x(i+2))-f(x(i)));
end;
b(n)=(3/h)*(f(x(n+2))-f(x(n)))-c1;
//
//resolution
//
p=A\b;
//
p=[c0;p;c1];
//
t=[];sp=[];
//
// reconstruction du profil
//
function s=reconst(t)
i=0;
while t>x(i+1)
i=i+1
end
if (t==0) then
   s=f(x(1))
elseif (t==1)
   s=f(x(n+2));
else
  xi=t-x(i);
  s=f(x(i))*(2*xi^3-3*h*xi^2+h^3)/h^3+f(x(i+1))*(-2*xi^3+3*h*xi^2)/h^3+...
    p(i)*(xi^3-2*h*xi^2+h^2*xi)/h^2+p(i+1)*(xi^3-h*xi^2)/h^2;	
end
//
endfunction
//
for i=1:50*n
  t(i)=(i-1)/(50*n-1);
  sp(i)=reconst(t(i));
end
//
// trace
//
clf()
plot2d([t,t],[f(t),sp],[1,2])
//
write(%io(2),'temps de calcul:');disp(timer())
//

