-- Solve for coefficients of a polynomial that goes through n data points.

include "::incl:solver"

data = read(xydata)
plot data[1:n]

x[i]=data[i,1] dim[n] -- separate x and y
y[i]=data[i,2] dim[n]

f(x)[j] = x^(j-1)  -- polynomial

-- construct matrix using x data points
A[i,j] = f(x[i])[j] dim[n,n]

-- solve for coeffs that give y points
c:=solve(A,y):
c:{0.797,3.863,-0.715,0.057,-0.002,0.000}

n=6
plot sum((c*f(X))[i],i,1,n)
