function [] = tri_quad();

% define nodal values
z = [
    0 0
    1 0
    0 1
    0.5 0.5
    0 0.5
    0.5 0
    ];

% generate 2d mesh for plotting
x = 0:0.01:1;
y = x;
[x,y] = meshgrid(x,y);

figure(1), clf;

% generate matrix;
A = zeros(6,6);
A(:,1) = ones(6,1);
A(:,2) = z(:,1);
A(:,3) = z(:,2);
A(:,4) = z(:,1).^2;
A(:,5) = z(:,1).*z(:,2);
A(:,6) = z(:,2).^2;

% six basis functions to compute and plot
for k = 1:6
    e = zeros(6,1);
    e(k) = 1;
    
    c = A\e;
    
    phi{k} = c(1) + ...
             c(2)*x + ...
             c(3)*y + ...
             c(4)*x.^2 + ...
             c(5)*x.*y + ...
             c(6)*y.^2;
    
    phi{k}(y > 1-x) = NaN;
    
    figure(1);
    subplot(3,2,k);
    surf(x,y,phi{k});
    shading interp
    set(gca,'FontSize',24);
    xlabel('x','FontSize',24);
    ylabel('y','FontSize',24);
    title(['\phi_',num2str(k)],'FontSize',24);
end

error('break');

% use these basis functions to interpolate exp(xy);
original = exp(x.*y);
original(y > 1-x) = NaN;
figure(2); clf;
subplot(1,2,1);
surf(x,y,original);
shading interp;
xlabel('x');
ylabel('y');
title('original');

% interpolated value
interp = zeros(size(x));
for k = 1:6
    uz = exp(z(k,1)*z(k,2));
    interp = interp + uz*phi{k};
end
figure(2)
subplot(1,2,2);
surf(x,y,original);
shading interp;
xlabel('x');
ylabel('y');
title('quadratic interpolation');