function cost = stamps(year)
% STAMPS  Fit high degree polynomial in Newton form through
%         stamp cost data.  Evaluates it at the input year.

stampdata  % fills t,y with data

ts = t - 1885;  % = [0 32 34 47 ... 106 110 114 116]

% build Newton matrix to find c_j in
% p(x) = c_1 + c_2 (x-0) + c_3 (x-0)(x-32) + c_4 (x-0)(x-32)(x-34)
%         + ... + c_18 (x-0)(x-32)...(x-114)
A = zeros(18,18);
A(1,1) = 1;
for i=2:18   % i = row index
  prodx = 1;
  for j=1:i  % j = col index
    A(i,j) = prodx;
    prodx = prodx * (ts(i)-ts(j));
  end
end
% to view checkable small chunk:  A(1:3,1:3)

c = A \ y;  % find coeffs: automatically recognizes triangular
            %   system and solves by forward substitution

% eval at input
x = year - 1885;
cost = c(1);
prodx = 1;
for j=2:18
  prodx = prodx * (x - ts(j-1));
  cost = cost + c(j) * prodx;
end
