You have to increase the number of values in x, when calling the function.
x=[0 8 16 25 37];
y=[5 17 12 21 19];
disp("Los datos del largo de la pieza son:")
x
disp('Los datos del ancho de la pieza son:')
y
n=length(x)-1;
xo=19;
disp("El orden de este polinomio es de: "+n)
%Diferencias divididas%
m1=(y(2)-y(1))/(x(2)-x(1));
m2=(y(3)-y(2))/(x(3)-x(2));
m3=(y(4)-y(3))/(x(4)-x(3));
m4=(y(5)-y(4))/(x(5)-x(4));
disp('')
%%%%%%%%%%%%%%%%%%%%%%%
b0=y(1);
b1=m1;
b2=(m2-b1)/(x(3)-x(1));
b3=(((m3-m2)/(x(4)-x(2))-b2)/(x(4)-x(1)));
b4=(((m4-m3)/(x(5)-x(3))-((m3-m2)/(x(4)-x(2))))/(x(5)-x(2))-(((m3-m2)/(x(4)-x(2)))-((m2-m1)/(x(3)-x(1))))/(x(4)-x(1)))/(x(5)-x(1));
fn= b0+b1*(xo-x(1))+b2*(xo-x(1))*(xo-x(2))+b3*(xo-x(1))*(xo-x(2))*(xo-x(3))+b4*(xo-x(1))*(xo-x(2))*(xo-x(3))*(xo-x(4));
disp("El resultado del polinomio es: "+fn);
plot(x,y,'p')
x = linspace(min(x),max(x)) ;
f4=@(x) 5+1.5*(x-0)-0.13281*(x-0).*(x-8)+0.009136*(x-0).*(x-8).*(x-16)-0.00038778*(x-0).*(x-8).*(x-16).*(x-25);
g=min(x)-3:1:max(x)+3;
hold
plot(x,f4(x))
grid
