clc;clear
format short e
w = [ 1 2 4 8 12 16 20 24 26 ];
T = [ 0 10 25 40 ];
rho = [ 1.00747, 1.00707, 1.00409, 0.99908 ; ...
        1.01509, 1.01442, 1.01112, 1.00593 ; ...
        1.03038, 1.02920, 1.02530, 1.01977 ; ...
        1.06121, 1.05907, 1.05412, 1.04798 ; ...
        1.09244, 1.08946, 1.08365, 1.07699 ; ...
        1.12419, 1.12056, 1.11401, 1.10688 ; ...
        1.15663, 1.15254, 1.14533, 1.13774 ; ...
        1.18999, 1.18577, 1.17776, 1.16971 ; ...
        1.20709, 1.20254, 1.19443, 1.18614 ] ;
n1 = length(w);
n2 = length(T);
for i = 1:n1 % create stacked vectors of independent variables
    for j = 1:n2
        wn((i-1)*n2+j) = w(i);
        Tn((i-1)*n2+j) = T(j);
    end
end
wn = wn';
Tn = Tn';
X = [ ones(n1*n2,1) wn Tn wn.^2 Tn.^2 wn.*Tn];  % create the X matrix
y = reshape(rho',[n1*n2,1]);  % create the y vector
A = X'*X;
rs = X'*y;
b = A\rs  % model coefficients
rhop = X*b;  % predicted densities
resid = y - rhop;  % residuals
plot(y,rhop,'ok')
grid
hold on
plot(rhop,rhop,'k--')
xlabel('Density Data - g/cm3')
ylabel('Predicted Density - g/cm3')
figure
plot(rhop,resid,'k.')
grid
xlabel('Predicted Density')
ylabel('Residual')
title('Residuals vs. Fits')
