function[] = oned()

    clear all;
    close all;

    % Create your geometry; not too difficult in 1D.
    L = 10;
    A = pi*(0.1^2);

    % Create the boundary tractions.  Since this is 1D, we just need the
    % tractions at the left and right ends.
    Fbar = [0; 500];

    % For Dirichelt BCs, we can't just have a list of values; we need to know which
    % nodes (or positions) have Dirichlet BCs, and then for each node (or position),
    % we need to know the value of the Dirichlet BC.
    %
    % Since we don't know which node corresponds to what physical location,
    % the positions in the DBC structures correspond to the left and right of the
    % domain.
    DBCNodeCode = [1; 0];
    DBCNodeValue = [0; 0];

    % Number of elements
    % You might find it interesting that I ran this little code with
    % 20,000+ elements (on a machine with enough memory), and it
    % returned a very nice answer.
    nel = 10;

    meshtype = 1; % 0 = typical
                  % 1 = random

    % Create the nodes, i.e., create the x positions of the nodes.
    xhat = Buildxhat(meshtype,L,nel);

    % Create the LM matrix; see the function below.
    LM = BuildLM(meshtype,L,nel,xhat);

    % Now, calculate the element stiffness matrix and force vector.
    K = zeros(???);
    F = zeros(???);
    % Loop over the elements.
    for iel = 1:nel
        % Assign the left and right endpoints of the element to a and b.
        a = xhat(LM(1,iel));
        b = xhat(LM(2,iel));

        % Zero out the element stiffness and force matrices.
        Ke = zeros(2,2);
        Fe = zeros(2,1);

        % I use Simpson's rule; look it up on Wikipedia.
        % Because of this, I have to do three evaluations, and then
        % a few ancillary calculations.
        %
        % If you were using, say, Gaussian quadrature, with the
        % integration points and weights in a nice data structure,
        % you might want a loop here.
        %
        % Note how the form of the Ke and Fe you calculated on paper shows up
        % here.
        [Ne,dNedx] = CalcInterpFunctions(a,b,a);
        E = CalcE(a);
        body = CalcBodyForce(a);
        % Note how I create Ke and Fe to be the correct shapes.
        % Ah, the beauty of matlab . . .
        Ke = Ke + A*E*dNedx'*dNedx;
        Fe = Fe + Ne'*body;

        [Ne,dNedx] = CalcInterpFunctions(a,b,(a+b)/2);
        E = CalcE((a+b)/2);
        body = CalcBodyForce((a+b)/2);
        Ke = Ke + 4*A*E*dNedx'*dNedx;
        Fe = Fe + 4*Ne'*body;

        [Ne,dNedx] = CalcInterpFunctions(a,b,b);
        E = CalcE(b);
        body = CalcBodyForce(b);
        Ke = Ke + A*E*dNedx'*dNedx;
        Fe = Fe + Ne'*body;

        Ke = (b-a)/6*Ke;
        Fe = (b-a)/6*Fe;

        % The tractions are just evaluations at points;
        % no integration needed, but I do need some conditional logic to
        % ensure that the evaluation is done properly . . .
        if (a == 0)
            [Ne,dNedx] = CalcInterpFunctions(a,b,0);
            Fe = Fe + Ne'*Fbar(1);
        elseif (b == L)
            [Ne,dNedx] = CalcInterpFunctions(a,b,L);
            Fe = Fe + Ne'*Fbar(2);
        end

        %Ke
        %Fe

        % Map local components to global components properly via LM.
        for i = 1:2
            for j = 1:2
                K = SOME FUNCTION OF Ke
            end
            F = SOME FUNCTION OF Fe
        end
    end

    % Deal with the Dirichlet BCs.
    Kmod = K;
    Fmod = F;
    if (DBCNodeCode(1) == 1)
        for inode = 1:nel+1
            if (xhat(inode) == 0)
                DBCNodeNum = inode;
                break;
            end
        end
        Fmod = Fmod - DBCNodeValue(1)*K(:,DBCNodeNum);
        Kmod(DBCNodeNum,:) = 0;
        Kmod(:,DBCNodeNum) = 0;
        Kmod(DBCNodeNum,DBCNodeNum) = 1;
        Fmod(DBCNodeNum,1) = DBCNodeValue(1);
    end
    if (DBCNodeCode(2) == 1)
        for inode = 1:nel+1
            if (xhat(inode) == L)
                DBCNodeNum = inode;
                break;
            end
        end
        Fmod = Fmod - DBCNodeValue(2)*K(:,DBCNodeNum);
        Kmod(DBCNodeNum,:) = 0;
        Kmod(:,DBCNodeNum) = 0;
        Kmod(DBCNodeNum,DBCNodeNum) = 1;
        Fmod(DBCNodeNum,1) = DBCNodeValue(2);
    end

    % Solve the system.
    uhat = Kmod\Fmod;

    % exact solution
    Ebar = 7e10;
    b0 = 500;
    b1 = 10;
    x = linspace(0,L);
    term1 = -log(1+x) + (1+x);
    term2 = (x.^2).*log(1+x) ...
           -((1+x).^2).*log(1+x) ...
           +((1+x).^2)/2 ...
           +2*(1+x).*log(1+x) ...
           -2*(1+x);
    term3 = log(1+x);
    c1 = 500 + b0*L + b1/2*L^2;
    c2 = 1/(A*Ebar)*(b0 - 3/4*b1);
    uexact = 1/(A*Ebar)*(-b0*term1 - b1/2*term2 + c1*term3) + c2;

    % Plot data for elements.
    for iel = 1:nel
        node1 = LM(1,iel);
        xhatplot(iel) = xhat(node1);
        uhatplot(iel) = uhat(node1);
    end
    node2 = LM(2,nel);
    xhatplot(nel+1) = xhat(node2);
    uhatplot(nel+1) = uhat(node2);

    % Plot it.
    plot(x,uexact,'b-');
    hold on;
    plot(xhatplot,uhatplot,'ro-');
    legend('exact','FE', ...
           'Location','best');

end

% Build the vector of nodal positions.
function[xhat] = Buildxhat(type,L,nel)
    if (type == 0)
        xhat = linspace(0,L,nel+1);
    else
        used = zeros(nel+1,1);
        nums = round((nel+1).*rand((nel+1)*100,1));
        i = 1;
        j = 0;
        h = L/nel;
        while (1)
            nodenum = nums(i);
            if (nodenum == 0)
                nodenum = 1;
            end
            i = i+1;
            if (used(nodenum)== 1)
                continue;
            end
            used(nodenum) = 1;
            xhat(nodenum) = j*h;
            j = j+1;
            if (j == nel+1)
                break;
            end
        end
    end
end

% Build the LM matrix.
function[LM] = BuildLM(type,L,nel,xhat)

    if (type == 0)
        for iel = 1:nel
            LM(1,iel) = iel;
            LM(2,iel) = iel+1;
        end
    else
        % The assumption here is that the elements are monotonically numbered,
        % with the left end of element 1 being at x=0, and the right end of
        % element nel being at x=L.
        h = L/nel;
        for iel = 1:nel
            a = h*(iel-1);
            b = h*iel;

            [STUFF DELETED]

        end
 
    end

end

% This returns the Young's modulus at any point x.
function[E] = CalcE(x)
    E = 7e10*(1+x);
end

% This returns the body force at any point x.
function[body] = CalcBodyForce(x)
    body = 500 + 10*x;
end

% This function calculates the interpolation functions at any point x,
% for the element with endpoints a < b.
% Note that they get stacked nicely into the array N, and likewise
% for the derivatives.
function[Ne,dNedx] = CalcInterpFunctions(a,b,x)
    N1e = (x-b)/(a-b);
    N2e = (x-a)/(b-a);
    Ne = [N1e N2e];

    dN1edx = 1/(a-b);
    dN2edx = 1/(b-a);
    dNedx = [dN1edx dN2edx];
end


