    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Deal with the Dirichlet BCs.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    Kmod = K;
    Fmod = F;
    for ibc = 1:4
        if (BCCode(ibc) == 2)
            continue;
        end
        % Figure out which nodes apply to the given edge.
        % There is no easy way to do this; check cases.
        if (ibc == 1)
            BCNodes = linspace(1,nel(2)*(nel(1)+1)+1,nel(2)+1);
        elseif (ibc == 2)
            BCNodes = nel(1) + linspace(1,nel(2)*(nel(1)+1)+1,nel(2)+1);
        elseif (ibc == 3)
            BCNodes = linspace(1,nel(1)+1,nel(1)+1);
        elseif (ibc == 4)
            BCNodes = nel(2)*(nel(1)+1) + linspace(1,nel(1)+1,nel(1)+1);
        end

        % loop over the nodes, calculate the BC at each node, and apply
        nBCNodes = max(size(BCNodes));
        for inode = 1:nBCNodes
            NodeNum = BCNodes(inode);
            xe = xhat(NodeNum,1);
            ye = xhat(NodeNum,2);
            BCVal = eval(char(DBCstr(ibc)));
            Fmod = Fmod - BCVal*Kmod(:,NodeNum);
            Kmod(NodeNum,:) = 0;
            Kmod(:,NodeNum) = 0;
            Kmod(NodeNum,NodeNum) = 1;
            Fmod(NodeNum,1) = BCVal;
        end
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Create plot data for elements.
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    for jnp = 1:(nel(2)+1)
        for inp = 1:(nel(1)+1)
            nodenum = (jnp-1)*(nel(1)+1) + inp;
            xplot(inp,jnp) = xhat(nodenum,1);
            yplot(inp,jnp) = xhat(nodenum,2);
            uplot(inp,jnp) = uhat(nodenum);
        end
    end

    % Plot it.
    contourf(xplot,yplot,uplot,100,'EdgeColor','none');
    colorbar('Location','Eastoutside');
    title(sprintf('nel = %d, l2 error = %e',nel(1)*nel(2),errorval));
    xlabel('x');
    ylabel('y');
    hold on;

    %print -dpdf me180_lab08.pdf;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function[Ne,dNedx,J] = CalcInterpFunctions(xi,eta,xhate,nen)
    % where:
    %   xi: the first isoparametric coordinate
    %   eta: the second isoparametric coordinate
    %   xhate: a vector of the nodal coordinates over a single element
    %   nen: the number of element nodes

    xi1 = -1;
    xi2 = 1;
    eta1 = -1;
    eta2 = 1;

    N1e = ???
    N2e = ???
    N3e = ???
    N4e = ???
    Ne = [N1e N2e N3e N4e];

    % Calc dNdxi
    dN1edxi = ???
    dN2edxi = ???
    dN3edxi = ???
    dN4edxi = ???
    dN1edeta = ???
    dN2edeta = ???
    dN3edeta = ???
    dN4edeta = ???
    dNedxi = [dN1edxi dN2edxi dN3edxi dN4edxi;
              dN1edeta dN2edeta dN3edeta dN4edeta];

    % Calc dNdx
    dxdxi = ???
    dydxi = ???
    dxdeta = ???
    dydeta = ???
    Jmat = [dxdxi dydxi;
            dxdeta dydeta];
    for ien = 1:nen
        dNiedxi = dNedxi(:,ien);
        dNiedx = Jmat\dNiedxi;
        dNedx(:,ien) = dNiedx;
    end

    % Calc J
    J = det(Jmat);

end


