function S = Spline3(a,b,x,y,N)
// Interpolation par splines cubiques 
// Voir pour les notations sur la feuille de TP jointe
    n  = length(x)
    dx = diff(x); dy = diff(y); dyodx = diff(y)./diff(x) 
    t = a + (b-a)*[0:1/N:1] 
// Construction de la matrice du systeme lineaire 
// dont sont solution les zk pour k =2 : n-1 
    M = zeros(n-2,n-2)
    M(1,1) = (dx(1) + dx(2))/3; M(1,2) = dx(2)/6 
    M(n-2,n-3) = dx(n-2)/6; M(n-2,n-2) = (dx(n-2) + dx(n-1))/6
    if n>4 then
        for k=2:n-3
            M(k,k-1) = dx(k)/6 
            M(k,k)   = (dx(k)+dx(k+1))/3
            M(k,k+1) = dx(k+1)/6
        end
    end
// Construction du second membre du systeme lineaire
// dont sont solutions les zk pour k=2 : n-1
   Y = []
   for k=1:n-2 
       Y = [Y; dy(k+1)/dx(k+1)-dy(k)/dx(k)]
   end       
// Calcul des zj pour j=1:n (en colonne)
   Z = inv(M)*Y; z = [0;Z;0]; dz = diff(z)  
// Calcul de la liste A des coefficients aj pour j=1:n-1 (en colonne)
   A = y(1:n-1)' - z(1:n-1).*(dx'.^2)/6    
// Calcul de la liste B des coefficients bj pour j=1:n-1 (en colonne)
   B = dyodx' - dz.*dx'/6  
// Calcul de la fonction spline sur le segment [x1,xn]
S  = zeros(t)
   for k= 1:n-1 
       TRONC = (sign((t-x(k)).*(x(k+1)-t))+ 1)/2 
       S    = S + (z(k+1)/dx(k))*TRONC.*((t-x(k)).^3)/6
       S    = S + (z(k)/dx(k))*TRONC.*((x(k+1)-t).^3)/6 
       S    = S + B(k)*TRONC.*(t-x(k)) + A(k)*TRONC
   end 
//Calcul de la fonction spline sur [a,x1]
       TRONC0 = (sign(x(1)-t)+ 1)/2
       S0     = TRONC0.*(y(1)+B(1)*(t-x(1)))
//Calcul de la fonction spline sur [xn,b] 
       TRONCn = (sign(t-x(n))+ 1)/2
       Sn     = TRONCn.*(y(n)+B(n-1)*(t-x(n)))
// Finalisation 
// Affichage graphique 
       plot(t,S+S0+Sn)
       set(gca(),"auto_clear","off")
       plot(x,y,'or')
       set(gca(),"auto_clear","on")                                                          
endfunction
