function [ret, msg] = interpNewtonCoefs(x, y)
  n = size(x, 1)
  if n <> size(y, 1) then
    [ret, msg] = return(0, "n <> size(y, 1)")
  end
  a = y
  for k = 2 : n
    for j = n : -1 : k
      a(j) = a(j) - a(j - 1)
      d = x(j) - x(j - k + 1)
      if d == 0 then
        [ret, msg] = return(0, "x(j) == x(j - k)")
      end
      a(j) = a(j) / d
    end
  end
  ret = zeros(n, 1)
  ret(1) = a(n)
  for j = n - 1 : -1 : 1
    for i = n - j + 1 : -1 : 2
      ret(i) = ret(i - 1) - ret(i) * x(j)
    end
    ret(1) = a(j) - ret(1) * x(j)
  end
  msg = "ok"
endfunction

function [ret, msg] = interpNewtonPoly(x, y, vname)
  [coeff, msg] = interpNewtonCoefs(x, y)
  if msg <> "ok" then
    [ret, msg] = return(0, msg)
  end
  ret = poly(coeff, vname, "coeff")
endfunction

x = [0.351;0.867;3.315;5.013;6.432]
y = [-0.572;-2.015;-3.342;-5.752;-6.911]
[interpPoly, msg] = interpNewtonPoly(x, y, "s")
 interpPoly = 

                                 2           3           4
   1.085347 -5.681919s +2.957907s  -0.650112s  +0.046258s