Scilab Program - divided difference formula - IProgramX

Program

 function[]=div(X,Y,x0)
      n=length(X)
      h=X(2)-X(1);
      for i=1:n
      Z(i)=Y(i);
      d(i,1)=Y(i);
  end
  for i=1:n-1
      for j=1:n-i
          d(i,j)=Z(j+1)-Z(j)/(X(j+i)-X(j));
          D(j,i+1)=d(i,j);
      end
      for k=1:n-i
      Z(k)=d(i,k);
  end
  end
  disp(D)
  p=Y(1);
  x=poly(0,'x')
  for i=1:n-1
       L=1;
       for j=1:i
           L=L*(x-X(j));
       end
       p=p+L*d(i,1);
   end
   disp(p,'newtons poly=')
  printf('value of poly @ %.4f=%.4f',x0,horner(p,x0))
  endfunction

Output

  exec('C:\Users\acer\div.sce', -1)

  X=[-1 1 4 6]
 X  =

  -1.   1.   4.   6.


  Y=[-4 8 -41 78]
 Y  =

  -4.   8.  -41.   78.


  div(X,Y,1.1)

   0.   10.        -45.666667   113.75714
   0.  -43.666667   107.23333   0.    
   0.   98.5        0.          0.    

 newtons poly=

                                    2            3
   506.69524 -103.75714x -500.69524x  +113.75714x
value of poly @ 1.1000=-61.8681