interpolacion polnomial

% interpolacion polinomial

x=[-1;0;1;2];
y=[-2;1;3;13];
n=length(x)-1;
b(:,n+1)=ones(size(x));
b(:,n)=x;

for k=n-1:-1:1
b(:,k)=b(:,k+1).*x;
end

b

a=inv(b)*y

%interpolacion

x0=[0.5;1.5];
y0=zeros(size(x0));

for k=1: n+1
y0=(y0+a(k)*x0.^(1+n-k)); k=1:n+1
end
y0

%grafica
xp=-1:0.01:2;
yp=zeros(size(xp));

for k=1:n+1
yp=(yp+a(k)*xp.^(1+n-k));

end

plot (xp,yp,x,y,'ob')
hold on
plot(x0,y0,'+r')