(* ::Package:: *)

  function a = Quintic(d,y)
% d is a column 4-vector with distinct components
% y is a column 4-vector
% a is a column 6-vector with the property that if
%
%    p(d) = a(1) + a(2)d + a(3)d^2 + a(4)d^3 + a(5)d^4 + a(6)d^5
%
% then  
%            p(d(1)) = y(1)
%            p(d(2)) = y(2)
%            p(d(3)) = y(3)
%            p(d(4)) = y(4)
%            p'(d(1)) = p'(d(4))
%            p "(d(1)) = p"(d(4))

% The 6 constraints are linear in a(1:6):

%        a(1) + d(1)*a(2) + d (1)^2*a(3) + d (1)^3*a(4) + d (1)^4*a(5) + d (1)^5*a(6) = y(1)
%        a(1) + d(2)*a(2) + d (2)^2*a(3) + d (2)^3*a(4) + d (2)^4*a(5) + d (2)^5*a(6) = y(2)
%        a(1) + d(3)*a(2) + d (3)^2*a(3) + d (3)^3*a(4) + d (3)^4*a(5) + d (3)^5*a(6) = y(3)
%        a(1) + d(4)*a(2) + d (4)^2*a(3) + d (4)^4*a(4) + d (4)^4*a(5) + d (4)^5*a(6) = y(4)
%        a(2) + 2*d(1)a(3) + 3*d (1)^2*a(4) + 4*d (1)^3*a(5) +5*d (1)^4 =
%                    a(2) + 2*d(4)a(3) + 3*d (4)^2*a(4) + 4*d (4)^3*a(5) = 5*d (4)^4 
%        2a(3) + 6d(1)*a(4) + 12d (1)^2*a(5) + 20*d (1)^3 = 
%               2a(3) + 6d(4)*a(4) + 12d (4)^2*a(5) + 20*d (4)^3

% Set up and solve the 6-by-6 linear system

a = ...
    [  1    d(1)     d (1)^2            d (1)^3            d (1)^4               d (1)^5;...
       1    d(2)     d (2)^2            d (2)^3            d (2)^4               d (2)^5;...
       1    d(3)     d (3)^2            d (3)^3            d (3)^4               d (3)^5;...
       1    d(4)     d (4)^2            d (4)^3            d (4)^4               d (4)^5;...
       0     0    2*(d(1)-d(4))  3*(d (1)^2-d (4)^2)   4*(d (1)^3-d (4)^3)    5*(d (1)^4 - d (4)^4);...
       0     0         0        6*(d(1)-d(4))       12*(d (1)^2-d (4)^2)  20*(d (1)^3-d (4)^3)   ] \ [y;0;0];



% If all errors are less than 1e-10: 5 points
% If the errors are larger than 1e-10 but still reasonable: 4 points
% If only the function values are correct but p' and p'' are not continuous at endpoints: 2 points
