function [M,polinomio]=Newton(xp,yp,xinterp)
  n = length(xp);
  
  M = zeros(n);
  
  M(:,1)=yp;
  for k=2:n
    for j=k:n
      M(j,k) = [M(j,k-1)-M(j-1,k-1)]/[xp(j)-xp(j-k+1)];
    endfor
  endfor
 % syms x
   polinomio=M(1,1);
   P=1;
   for k=1:n-1
     P = P*(xinterp-xp(k));
     polinomio = polinomio+P*M(k+1, k+1);
   endfor
  %polnew = expand(polnew);
  %pretty(polnew)
