Interpolacion spline segundo grado

function [A,D,sol]=spline2(x,y)
%Funcion basada en el metodo de Spline de 2do grado con los siguientes
%parametros:
%x:vector en x
%y: vector f(x)
% Y la funcion arroja los siguientes resultados:
%A: matriz cuadrada de coeficientes
%D: vector demismas dimensiones
%sol: Solucion Sistema por gauss
n=length(x);
m=(n-1)*3-1;
r=(n*2)-3;
A=zeros(m,m);
D=zeros(m,1);
l=3;p=1;
for i=1:2
for j=1:2
A(i,j)=x(i).^(-j+2);
end
end

for i=3:2:r
for j=i:i+1
p=j-p;for k=-2:0
A(j,((l+2)+k))=x(p).^(-k);
end
p=p-1;
end
l=l+3;
end
p=3;
j=3;l=1;
A(r+2,1)=1;
for i=3:4
A(r+2,i)=-((2*x(2)).^(3-(i-1)));
end
for i=n*2:m
for j=j:3:j+3
for k=j:j+1A(i,k)=l*((2*x(p)).^(j-(k-1)));
end
l=l*(-1);
end
p=p+1;
end
D(1,1)=y(1);
j=2;
fori=2:n-1
for j=j:j+1
D(j,1)=y(i);
end
j=j+1;
end
D(r+1,1)=y(n);
sol=gausselim(A,D);
plot(x,y,'or'), hold onf= @(x)(sol(1)*x+sol(2));
fplot(f,[x(1) x(2)]), hold on
j=3;
for i=3:3:m-2
f= @(x)(sol(i)*x^2+sol(i+1)*x+sol(i+2));fplot(f,[x(j-1) x(j)]), hold on
j=j+1;
end
xlim([x(1) x(n)])
title('Interpolación Splin Cuadrática')
hold off;