Coordinatematrix = 0;
Connectivitymatrix = 0;
Lengthmatrix = 0;
Lambdamatrix = 0;
Degreesoffreedommatrix = 0;
Localstiffnessmatrix3 = 0;
Localstiffnessmatrix4 = 0;
Ans1 = 0;
Ans2 = 0;
Ans3 = 0;
Ans4 = 0;
n_node = input('Input number of nodes');
n_element = input('Input number of elements');
n_dof_node = input('Input number of degrees of freedom per node(3 for frame/beam, 2 for truss)');
n_node_element = 2;
n_dof = n_dof_node*n_node;
reply = input('Please enter if your problem is a Truss or Frame/Beam T/F [T]: ','s');
if isempty(reply) ;
    reply = 'T';
end  ;
Areamatrix = zeros(n_element,1);
Elasticmodulusmatrix = zeros(n_element,1);
Restraintmatrix = zeros(n_dof,1);
Momentofinertiamatrix = zeros(n_element,1);
Forcematrix = zeros(n_dof,1);
Forcematrix2 = zeros(n_dof,1);
Displacementmatrix = zeros(n_dof,1); 
Displacementmatrix2 = zeros(n_dof,1); 
if reply == 'T';
    element = 1;
    while element <= n_element;
        Areamatrix (element,1)= input('Input area of element in (m) or (in)');
        Elasticmodulusmatrix (element,1)= input('Input modulus of elasticity of element in (kPa) or (ksi)');
        element = element + 1;
    end;
    node = 1;
    while node <= n_node;
        Coordinatematrix (node,1) = input('Input global x coordinate of node in (m) or (in)');
        Coordinatematrix (node,2) = input('Input global y coordinate of node in (m) or (in)');
        Degreesoffreedommatrix (node,1) = input('Input # of x DOF of node');
        Degreesoffreedommatrix (node,2) = input('Input # of y DOF of node');
        Ans = input('Is there a hinge at this node? [ 1 / 0 ]');
        if Ans == 1 ;
            Restraintmatrix (Degreesoffreedommatrix(node,1),1) = Degreesoffreedommatrix(node,1);
            Restraintmatrix (Degreesoffreedommatrix(node,2), 1) = Degreesoffreedommatrix(node,2);
        else;
            Ans2 = input('Is there a roller at this node? [ 1 / 0 ]');
            if Ans2 == 1 ;
                Ans3 = input('Input the degree of freedom that is restrained');
                if Ans3 == Degreesoffreedommatrix(node,1);
                    Restraintmatrix(Degreesoffreedommatrix(node,1),1) = Degreesoffreedommatrix(node,1) ;
                    Displacementmatrix2(Degreesoffreedommatrix(node,2),1) = 1;
                elseif Ans3 == Degreesoffreedommatrix(node,2);
                    Restraintmatrix(Degreesoffreedommatrix(node,2),1) = Degreesoffreedommatrix(node,2)     ;
                    Displacementmatrix2(Degreesoffreedommatrix(node,1),1) = 1;
                end;
            else;
                Displacementmatrix2(Degreesoffreedommatrix(node,1),1) = 1;
                Displacementmatrix2(Degreesoffreedommatrix(node,2),1) = 1;
            end;
        end;
        Ans = input('Is there a force at this node? [ 1 / 0 ]');
        while Ans == 1;
            Ans2 = input('Imput magnitude of force. Up&right = +, down&left = -');
            Ans3 = input('Imput degree of freedom that force works along');
            if Ans3 == Degreesoffreedommatrix(node,1);
                if Restraintmatrix(Degreesoffreedommatrix(node,1),1) == 0;
                    Forcematrix(Degreesoffreedommatrix(node,1),1) = Ans2;
                else;
                    Forcematrix2(Degreesoffreedommatrix(node,1),1) = Ans2;
                end;
            end;
            if Ans3 == Degreesoffreedommatrix(node,2);
                if Restraintmatrix(Degreesoffreedommatrix(node,2),1) == 0;
                    Forcematrix(Degreesoffreedommatrix(node,2),1) = Ans2;
                else;
                    Forcematrix2(Degreesoffreedommatrix(node,2),1) = Ans2;
                end;
            end;
            Ans = input('Is there another force at this node? [1 / 0]')       ;
        end       ;
        node = node + 1;
    end;
    element = 1;
    while element <= n_element;
        Connectivitymatrix (element, 1) = input('Input # of near node');
        Connectivitymatrix (element, 2) = element;
        Connectivitymatrix (element, 3) = input('Input # of far node');
        element = element + 1;
    end;
    element = 1;
    Localstiffnessmatrix  = zeros(n_dof, n_dof);
    Localstiffnessmatrix2 = zeros(n_dof, n_dof);
    while element <= n_element;
        nodenear = Connectivitymatrix (element,1);
        nodefar = Connectivitymatrix (element,3);
        nearx = Coordinatematrix (nodenear, 1);
        neary = Coordinatematrix (nodenear, 2);
        farx = Coordinatematrix (nodefar, 1);
        fary = Coordinatematrix (nodefar, 2);
        Lengthmatrix(element,1) = sqrt ((farx - nearx)^2+(fary-neary)^2);
        Lambdamatrix(element, 1) = (farx - nearx)/Lengthmatrix(element,1);
        Lambdamatrix(element, 2) = (fary - neary)/Lengthmatrix(element,1);
        Lambdamatrix(element, 3) = Lambdamatrix(element,1)^2;
        Lambdamatrix(element, 4) = Lambdamatrix(element,2)^2;
        Lambdamatrix(element, 5) = Lambdamatrix(element,1)*Lambdamatrix(element,2) ;
        Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,1), Degreesoffreedommatrix (nodenear,1)) = Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,1), Degreesoffreedommatrix (nodenear,1)) + Lambdamatrix(element,3)*Areamatrix (element,1)*Elasticmodulusmatrix (element,1)/Lengthmatrix (element, 1);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,1), Degreesoffreedommatrix (nodenear,2)) = Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,1), Degreesoffreedommatrix (nodenear,2)) + Lambdamatrix(element,5)*Areamatrix (element,1)*Elasticmodulusmatrix (element,1)/Lengthmatrix (element, 1);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,1), Degreesoffreedommatrix (nodefar,1)) = Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,1), Degreesoffreedommatrix (nodefar,1)) - Lambdamatrix(element,3)*Areamatrix (element,1)*Elasticmodulusmatrix (element,1)/Lengthmatrix (element, 1);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,1), Degreesoffreedommatrix (nodefar,2)) = Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,1), Degreesoffreedommatrix (nodefar,2)) - Lambdamatrix(element,5)*Areamatrix (element,1)*Elasticmodulusmatrix (element,1)/Lengthmatrix (element, 1);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,2), Degreesoffreedommatrix (nodenear,1)) = Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,2), Degreesoffreedommatrix (nodenear,1)) + Lambdamatrix(element,5)*Areamatrix (element,1)*Elasticmodulusmatrix (element,1)/Lengthmatrix (element, 1);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,2), Degreesoffreedommatrix (nodenear,2)) = Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,2), Degreesoffreedommatrix (nodenear,2)) + Lambdamatrix(element,4)*Areamatrix (element,1)*Elasticmodulusmatrix (element,1)/Lengthmatrix (element, 1);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,2), Degreesoffreedommatrix (nodefar,1)) =  Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,2), Degreesoffreedommatrix (nodefar,1)) - Lambdamatrix(element,5)*Areamatrix (element,1)*Elasticmodulusmatrix (element,1)/Lengthmatrix (element, 1);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,2), Degreesoffreedommatrix (nodefar,2)) = Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,2), Degreesoffreedommatrix (nodefar,2)) - Lambdamatrix(element,4)*Areamatrix (element,1)*Elasticmodulusmatrix (element,1)/Lengthmatrix (element, 1);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodefar,1), Degreesoffreedommatrix (nodenear,1)) = Localstiffnessmatrix(Degreesoffreedommatrix (nodefar,1), Degreesoffreedommatrix (nodenear,1)) - Lambdamatrix(element,3)*Areamatrix (element,1)*Elasticmodulusmatrix (element,1)/Lengthmatrix (element, 1);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodefar,1), Degreesoffreedommatrix (nodenear,2)) = Localstiffnessmatrix(Degreesoffreedommatrix (nodefar,1), Degreesoffreedommatrix (nodenear,2)) - Lambdamatrix(element,5)*Areamatrix (element,1)*Elasticmodulusmatrix (element,1)/Lengthmatrix (element, 1);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodefar,1), Degreesoffreedommatrix (nodefar,1)) = Localstiffnessmatrix( Degreesoffreedommatrix (nodefar,1), Degreesoffreedommatrix (nodefar,1)) + Lambdamatrix(element,3)*Areamatrix (element,1)*Elasticmodulusmatrix (element,1)/Lengthmatrix (element, 1);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodefar,1), Degreesoffreedommatrix (nodefar,2)) = Localstiffnessmatrix( Degreesoffreedommatrix (nodefar,1), Degreesoffreedommatrix (nodefar,2)) + Lambdamatrix(element,5)*Areamatrix (element,1)*Elasticmodulusmatrix (element,1)/Lengthmatrix (element, 1);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodefar,2), Degreesoffreedommatrix (nodenear,1)) = Localstiffnessmatrix( Degreesoffreedommatrix (nodefar,2), Degreesoffreedommatrix (nodenear,1)) - Lambdamatrix(element,5)*Areamatrix (element,1)*Elasticmodulusmatrix (element,1)/Lengthmatrix (element, 1);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodefar,2), Degreesoffreedommatrix (nodenear,2)) = Localstiffnessmatrix( Degreesoffreedommatrix (nodefar,2), Degreesoffreedommatrix (nodenear,2)) - Lambdamatrix(element,4)*Areamatrix (element,1)*Elasticmodulusmatrix (element,1)/Lengthmatrix (element, 1);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodefar,2), Degreesoffreedommatrix (nodefar,1)) =  Localstiffnessmatrix( Degreesoffreedommatrix (nodefar,2), Degreesoffreedommatrix (nodefar,1)) + Lambdamatrix(element,5)*Areamatrix (element,1)*Elasticmodulusmatrix (element,1)/Lengthmatrix (element, 1);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodefar,2), Degreesoffreedommatrix (nodefar,2)) = Localstiffnessmatrix( Degreesoffreedommatrix (nodefar,2), Degreesoffreedommatrix (nodefar,2)) + Lambdamatrix(element,4)*Areamatrix (element,1)*Elasticmodulusmatrix (element,1)/Lengthmatrix (element, 1);
        element = element + 1;
    end;
    counter = 1;
    Localstiffnessmatrix2 = Localstiffnessmatrix;
    while counter <= n_dof;
        if Displacementmatrix2(counter,1) == 0;
            Localstiffnessmatrix2(:,counter) = 0;
            Localstiffnessmatrix2(counter,:) = 0;
            Localstiffnessmatrix2(counter,counter) = 1;
        end;
        counter = counter + 1;
    end;
    Displacementmatrix2 = Localstiffnessmatrix2\Forcematrix;
    Forcematrix3 = Localstiffnessmatrix*Displacementmatrix2;
    Forcematrix3 = Forcematrix3 - Forcematrix2;
elseif reply == 'F';
    element = 1;
    while element <= n_element;
        Areamatrix (element,1)= input('Input area of element in (m) or (in)');
        Elasticmodulusmatrix (element,1)= input('Input modulus of elasticity of element in (kPa) or (ksi)');
        Momentofinertiamatrix (element,1)= input('Input moment of inertia of element in (m) or (in)');
        element = element + 1;
    end;
    node = 1;
    while node <= n_node;
        Coordinatematrix (node,1) = input('Imput global x coordinate of node in (m) or (in)');
        Coordinatematrix (node,2) = input('Imput global y coordinate of node in (m) or (in)');
        Degreesoffreedommatrix (node,1) = input('Imput # of x DOF of node');
        Degreesoffreedommatrix (node,2) = input('Imput # of y DOF of node');
        Degreesoffreedommatrix (node,3) = input('Imput # of Mz (i.e. rotation) DOF of node');
        Ans = input('Is it fixed at this node? [ 1 / 0 ]');
        if Ans == 1;
            Restraintmatrix(Degreesoffreedommatrix(node,1),1) = Degreesoffreedommatrix(node,1);
            Restraintmatrix(Degreesoffreedommatrix(node,2),1) = Degreesoffreedommatrix(node,2);
            Restraintmatrix(Degreesoffreedommatrix(node,3),1) = Degreesoffreedommatrix(node,3);
            Displacementmatrix2(Degreesoffreedommatrix(node,1),1) = 1;
            Displacementmatrix2(Degreesoffreedommatrix(node,2),1) = 1;
            Displacementmatrix2(Degreesoffreedommatrix(node,3,1)) = 1;
        elseif Ans == 0    ;
            Ans2 = input('Is there a hinge at this node? [ 1 / 0 ]');
            if Ans2 == 1 ;
                Restraintmatrix (Degreesoffreedommatrix(node,1),1) = Degreesoffreedommatrix(node,1);
                Restraintmatrix (Degreesoffreedommatrix(node,2), 1) = Degreesoffreedommatrix(node,2);
                Displacementmatrix2(Degreesoffreedommatrix(node,1),1) = 1;
                Displacementmatrix2(Degreesoffreedommatrix(node,2),1) = 1 ;
            elseif Ans2 == 0;
                Ans4 = input('Is there a roller at this node? [ 1 / 0 ]');
                if Ans4 == 1 ;
                    Ans3 = input('Input the degree of freedom that is restrained');
                    if Ans3 == Degreesoffreedommatrix(node,1);
                        Restraintmatrix(Degreesoffreedommatrix(node,1),1) = Degreesoffreedommatrix(node,1) ;
                        Displacementmatrix2(Degreesoffreedommatrix(node,1),1) = 1;
                    elseif Ans3 == Degreesoffreedommatrix(node,2);
                        Restraintmatrix(Degreesoffreedommatrix(node,2),1) = Degreesoffreedommatrix(node,2)      ;
                        Displacementmatrix2(Degreesoffreedommatrix(node,2),1) = 1;
                    end;
                end;
            end;
        end;
        Ans = input('Is there a force/moment at this node? [ 1 / 0 ]');
        while Ans == 1;
            Ans2 = input('Imput magnitude of force. Up&right = +, down&left = -');
            Ans3 = input('Imput degree of freedom that force works along');
            if Ans3 == Degreesoffreedommatrix(node,1);
                if Restraintmatrix(Degreesoffreedommatrix(node,1),1) == 0;
                    Forcematrix(Degreesoffreedommatrix(node,1),1) = Ans2;
                else ;
                    Forcematrix2(Degreesoffreedommatrix(node,1),1) = Ans2;
                end;
            end;
            if Ans3 == Degreesoffreedommatrix(node,2);
                if Restraintmatrix(Degreesoffreedommatrix(node,2),1) == 0;
                    Forcematrix(Degreesoffreedommatrix(node,2),1) = Ans2;
                else;
                    Forcematrix2(Degreesoffreedommatrix(node,2),1) = Ans2;
                end;
            end;
            if Ans3 == Degreesoffreedommatrix(node,3);
                if Restraintmatrix(Degreesoffreedommatrix(node,3),1) == 0;
                    Forcematrix(Degreesoffreedommatrix(node,3),1) = Ans2;
                else;
                    Forcematrix2(Degreesoffreedommatrix(node,3),1) = Ans2;
                end;
            end;
            Ans = input('Is there another force/moment at this node?');
        end    ; 
        node = node + 1;
    end;
    element = 1;
    while element <= n_element;
        Connectivitymatrix (element, 1) = input('Imput # of near node');
        Connectivitymatrix (element, 2) = element;
        Connectivitymatrix (element, 3) = input('Imput # of far node');
        element = element + 1;
    end;
    element = 1;
    Localstiffnessmatrix = zeros (n_dof, n_dof);
    while element <= n_element;
        nodenear = Connectivitymatrix (element,1);
        nodefar = Connectivitymatrix (element,3);
        nearx = Coordinatematrix (nodenear, 1);
        neary = Coordinatematrix (nodenear, 2);
        farx = Coordinatematrix (nodefar, 1);
        fary = Coordinatematrix (nodefar, 2);
        Lengthmatrix(element,1) = sqrt ((farx - nearx)^2+(fary-neary)^2);
        Lambdamatrix(element, 1) = (farx - nearx)/Lengthmatrix(element,1);
        Lambdamatrix(element, 2) = (fary - neary)/Lengthmatrix(element,1);
        Lambdamatrix(element, 3) = Lambdamatrix(element,1)^2;
        Lambdamatrix(element, 4) = Lambdamatrix(element,2)^2;
        Lambdamatrix(element, 5) = Lambdamatrix(element,1)*Lambdamatrix(element,2);
        A = Areamatrix (element,1);
        E = Elasticmodulusmatrix (element,1);
        L = Lengthmatrix (element, 1);
        I = Momentofinertiamatrix (element,1);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,1), Degreesoffreedommatrix (nodenear,1)) = Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,1), Degreesoffreedommatrix (nodenear,1)) + (A*E*Lambdamatrix(element, 3)/L + 12*E*I*Lambdamatrix(element, 4)/L^3);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,1), Degreesoffreedommatrix (nodenear,2)) = Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,1), Degreesoffreedommatrix (nodenear,2)) + (A*E/L - 12*E*I/L^3)*Lambdamatrix(element, 5);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,1), Degreesoffreedommatrix (nodenear,3)) = Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,1), Degreesoffreedommatrix (nodenear,3)) - (6*E*I/L^2)*Lambdamatrix(element,2);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,1), Degreesoffreedommatrix (nodefar,1)) = Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,1), Degreesoffreedommatrix (nodefar,1)) - (A*E*Lambdamatrix(element, 3)/L + 12*E*I*Lambdamatrix(element, 4)/L^3);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,1), Degreesoffreedommatrix (nodefar,2)) = Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,1), Degreesoffreedommatrix (nodefar,2)) - (A*E/L - 12*E*I/L^3)*Lambdamatrix(element, 5);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,1), Degreesoffreedommatrix (nodefar,3)) = Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,1), Degreesoffreedommatrix (nodefar,3)) - (6*E*I/L^2)*Lambdamatrix(element,2);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,2), Degreesoffreedommatrix (nodenear,1)) = Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,2), Degreesoffreedommatrix (nodenear,1)) + (A*E/L - 12*E*I/L^3)*Lambdamatrix(element, 5);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,2), Degreesoffreedommatrix (nodenear,2)) = Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,2), Degreesoffreedommatrix (nodenear,2)) + (A*E*Lambdamatrix(element, 4)/L + 12*E*I*Lambdamatrix(element, 3)/L^3);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,2), Degreesoffreedommatrix (nodenear,3)) = Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,2), Degreesoffreedommatrix (nodenear,3)) + (6*E*I/L^2)*Lambdamatrix(element,1)    ;
        Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,2), Degreesoffreedommatrix (nodefar,1)) =  Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,2), Degreesoffreedommatrix (nodefar,1)) - (A*E/L - 12*E*I/L^3)*Lambdamatrix(element, 5);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,2), Degreesoffreedommatrix (nodefar,2)) = Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,2), Degreesoffreedommatrix (nodefar,2)) - (A*E*Lambdamatrix(element, 4)/L + 12*E*I*Lambdamatrix(element, 3)/L^3);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,2), Degreesoffreedommatrix (nodefar,3)) = Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,2), Degreesoffreedommatrix (nodefar,3)) + (6*E*I/L^2)*Lambdamatrix(element,1)    ;
        Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,3), Degreesoffreedommatrix (nodenear,1)) = Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,3), Degreesoffreedommatrix (nodenear,1)) - (6*E*I/L^2)*Lambdamatrix(element,2);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,3), Degreesoffreedommatrix (nodenear,2)) = Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,3), Degreesoffreedommatrix (nodenear,2)) + (6*E*I/L^2)*Lambdamatrix(element,1);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,3), Degreesoffreedommatrix (nodenear,3)) = Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,3), Degreesoffreedommatrix (nodenear,3)) + (4*E*I/L);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,3), Degreesoffreedommatrix (nodefar,1)) = Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,3), Degreesoffreedommatrix (nodefar,1)) + (6*E*I/L^2)*Lambdamatrix(element,2);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,3), Degreesoffreedommatrix (nodefar,2)) = Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,3), Degreesoffreedommatrix (nodefar,2)) - (6*E*I/L^2)*Lambdamatrix(element,1);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,3), Degreesoffreedommatrix (nodefar,3)) = Localstiffnessmatrix(Degreesoffreedommatrix (nodenear,3), Degreesoffreedommatrix (nodefar,3)) + (2*E*I/L);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodefar,1), Degreesoffreedommatrix (nodenear,1)) = Localstiffnessmatrix(Degreesoffreedommatrix (nodefar,1), Degreesoffreedommatrix (nodenear,1)) - (A*E*Lambdamatrix(element, 3)/L + 12*E*I*Lambdamatrix(element, 4)/L^3);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodefar,1), Degreesoffreedommatrix (nodenear,2)) = Localstiffnessmatrix(Degreesoffreedommatrix (nodefar,1), Degreesoffreedommatrix (nodenear,2)) - (A*E/L - 12*E*I/L^3)*Lambdamatrix(element, 5);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodefar,1), Degreesoffreedommatrix (nodenear,3)) = Localstiffnessmatrix(Degreesoffreedommatrix (nodefar,1), Degreesoffreedommatrix (nodenear,3)) + (6*E*I/L^2)*Lambdamatrix(element,2);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodefar,1), Degreesoffreedommatrix (nodefar,1)) = Localstiffnessmatrix( Degreesoffreedommatrix (nodefar,1), Degreesoffreedommatrix (nodefar,1)) + (A*E*Lambdamatrix(element, 3)/L + 12*E*I*Lambdamatrix(element, 4)/L^3);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodefar,1), Degreesoffreedommatrix (nodefar,2)) = Localstiffnessmatrix( Degreesoffreedommatrix (nodefar,1), Degreesoffreedommatrix (nodefar,2)) + (A*E/L - 12*E*I/L^3)*Lambdamatrix(element, 5);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodefar,1), Degreesoffreedommatrix (nodefar,3)) = Localstiffnessmatrix( Degreesoffreedommatrix (nodefar,1), Degreesoffreedommatrix (nodefar,3)) + (6*E*I/L^2)*Lambdamatrix(element,2);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodefar,2), Degreesoffreedommatrix (nodenear,1)) = Localstiffnessmatrix( Degreesoffreedommatrix (nodefar,2), Degreesoffreedommatrix (nodenear,1)) - (A*E/L - 12*E*I/L^3)*Lambdamatrix(element, 5);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodefar,2), Degreesoffreedommatrix (nodenear,2)) = Localstiffnessmatrix( Degreesoffreedommatrix (nodefar,2), Degreesoffreedommatrix (nodenear,2)) - (A*E*Lambdamatrix(element, 4)/L + 12*E*I*Lambdamatrix(element, 3)/L^3);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodefar,2), Degreesoffreedommatrix (nodenear,3)) = Localstiffnessmatrix( Degreesoffreedommatrix (nodefar,2), Degreesoffreedommatrix (nodenear,3)) - (6*E*I/L^2)*Lambdamatrix(element,1);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodefar,2), Degreesoffreedommatrix (nodefar,1)) =  Localstiffnessmatrix( Degreesoffreedommatrix (nodefar,2), Degreesoffreedommatrix (nodefar,1)) + (A*E/L - 12*E*I/L^3)*Lambdamatrix(element, 5);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodefar,2), Degreesoffreedommatrix (nodefar,2)) = Localstiffnessmatrix( Degreesoffreedommatrix (nodefar,2), Degreesoffreedommatrix (nodefar,2)) + (A*E*Lambdamatrix(element, 4)/L + 12*E*I*Lambdamatrix(element, 3)/L^3);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodefar,2), Degreesoffreedommatrix (nodefar,3)) = Localstiffnessmatrix( Degreesoffreedommatrix (nodefar,2), Degreesoffreedommatrix (nodefar,3)) - (6*E*I/L^2)*Lambdamatrix(element,1);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodefar,3), Degreesoffreedommatrix (nodenear,1)) = Localstiffnessmatrix( Degreesoffreedommatrix (nodefar,3), Degreesoffreedommatrix (nodenear,1)) - (6*E*I/L^2)*Lambdamatrix(element,2);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodefar,3), Degreesoffreedommatrix (nodenear,2)) = Localstiffnessmatrix( Degreesoffreedommatrix (nodefar,3), Degreesoffreedommatrix (nodenear,2)) + (6*E*I/L^2)*Lambdamatrix(element,1);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodefar,3), Degreesoffreedommatrix (nodenear,3)) = Localstiffnessmatrix( Degreesoffreedommatrix (nodefar,3), Degreesoffreedommatrix (nodenear,3)) + (2*E*I/L);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodefar,3), Degreesoffreedommatrix (nodefar,1)) =  Localstiffnessmatrix( Degreesoffreedommatrix (nodefar,3), Degreesoffreedommatrix (nodefar,1)) + (6*E*I/L^2)*Lambdamatrix(element,2);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodefar,3), Degreesoffreedommatrix (nodefar,2)) = Localstiffnessmatrix( Degreesoffreedommatrix (nodefar,3), Degreesoffreedommatrix (nodefar,2)) - (6*E*I/L^2)*Lambdamatrix(element,1);
        Localstiffnessmatrix(Degreesoffreedommatrix (nodefar,3), Degreesoffreedommatrix (nodefar,3)) = Localstiffnessmatrix( Degreesoffreedommatrix (nodefar,3), Degreesoffreedommatrix (nodefar,3)) + (4*E*I/L);
        Ans = input('Is there a uniform load on this element? [ 1 / 0 ]');
        if Ans == 1;
            Ans2 = input('Enter the magnitude(weight/unitlength) of the load (+ up, - down)');
            if Restraintmatrix(Degreesoffreedommatrix(nodenear,2),1) == 0;
                Forcematrix(Degreesoffreedommatrix(nodenear,2),1) = Forcematrix(Degreesoffreedommatrix(nodenear,2),1) + Ans2*Lengthmatrix(element,1)/2;
            else;
                Forcematrix2(Degreesoffreedommatrix(nodenear,2),1) = Forcematrix2(Degreesoffreedommatrix(nodenear,2),1) + Ans2*Lengthmatrix(element,1)/2;
            end;
            if Restraintmatrix(Degreesoffreedommatrix(nodefar,2),1) == 0;
                Forcematrix(Degreesoffreedommatrix(nodefar,2),1) = Forcematrix(Degreesoffreedommatrix(nodefar,2),1) + Ans2*Lengthmatrix(element,1)/2;
            else;
                Forcematrix2(Degreesoffreedommatrix(nodefar,2),1) = Forcematrix2(Degreesoffreedommatrix(nodefar,2),1) + Ans2*Lengthmatrix(element,1)/2;
            end;
            if Restraintmatrix(Degreesoffreedommatrix(nodenear,3),1) == 0;
                Forcematrix(Degreesoffreedommatrix(nodenear,3),1) = Forcematrix(Degreesoffreedommatrix(nodenear,3),1) + Ans2*Lambdamatrix(element, 1)*(Lengthmatrix(element,1)^2)/12;
            else;
                Forcematrix2(Degreesoffreedommatrix(nodenear,3),1) = Forcematrix2(Degreesoffreedommatrix(nodenear,3),1) + Ans2*Lambdamatrix(element, 1)*(Lengthmatrix(element,1)^2)/12;
            end;
            if Restraintmatrix(Degreesoffreedommatrix(nodefar,3),1) == 0;
                Forcematrix(Degreesoffreedommatrix(nodefar,3),1) = Forcematrix(Degreesoffreedommatrix(nodefar,3),1) - Ans2*Lambdamatrix(element, 1)*(Lengthmatrix(element,1)^2)/12;
            else;
                Forcematrix2(Degreesoffreedommatrix(nodefar,3),1) = Forcematrix2(Degreesoffreedommatrix(nodefar,3),1) - Ans2*Lambdamatrix(element, 1)*(Lengthmatrix(element,1)^2)/12;
            end;
        end;
        Ans = input('Is there a point load on the middle of this element? [ 1 / 0 ]');
        if Ans == 1;
            Ans2 = input('Enter the magnitude(weight/unitlength) of the load (+ up, - down)');
            if Restraintmatrix(Degreesoffreedommatrix(nodenear,2),1) == 0;
                Forcematrix(Degreesoffreedommatrix(nodenear,2),1) = Forcematrix(Degreesoffreedommatrix(nodenear,2),1) + Ans2/2;
            else;
                Forcematrix2(Degreesoffreedommatrix(nodenear,2),1) = Forcematrix2(Degreesoffreedommatrix(nodenear,2),1) + Ans2/2;
            end;
            if Restraintmatrix(Degreesoffreedommatrix(nodefar,2),1) == 0;
                Forcematrix(Degreesoffreedommatrix(nodefar,2),1) = Forcematrix(Degreesoffreedommatrix(nodefar,2),1) + Ans2/2;
            else;
                Forcematrix2(Degreesoffreedommatrix(nodefar,2),1) = Forcematrix2(Degreesoffreedommatrix(nodefar,2),1) + Ans2/2;
            end;
            if Restraintmatrix(Degreesoffreedommatrix(nodenear,3),1) == 0;
                Forcematrix(Degreesoffreedommatrix(nodenear,3),1) = Forcematrix(Degreesoffreedommatrix(nodenear,3),1) + Ans2*Lambdamatrix(element, 1)*Lengthmatrix(element,1)/8;
            else;
                Forcematrix2(Degreesoffreedommatrix(nodenear,3),1) = Forcematrix2(Degreesoffreedommatrix(nodenear,3),1) + Ans2*Lambdamatrix(element, 1)*Lengthmatrix(element,1)/8;
            end;
            if Restraintmatrix(Degreesoffreedommatrix(nodefar,3),1) == 0;
                Forcematrix(Degreesoffreedommatrix(nodefar,3),1) = Forcematrix(Degreesoffreedommatrix(nodefar,3),1) - Ans2*Lambdamatrix(element, 1)*Lengthmatrix(element,1)/8;
            else;
                Forcematrix2(Degreesoffreedommatrix(nodefar,3),1) = Forcematrix2(Degreesoffreedommatrix(nodefar,3),1) - Ans2*Lambdamatrix(element, 1)*Lengthmatrix(element,1)/8;
            end;
        end;
        element = element + 1;
    end;
    counter = 1;
    Localstiffnessmatrix2 = Localstiffnessmatrix;
    while counter <= n_dof;
        if Displacementmatrix2(counter,1) == 1;
            Localstiffnessmatrix2(:,counter) = 0;
            Localstiffnessmatrix2(counter,:) = 0;
            Localstiffnessmatrix2(counter,counter) = 1;
        end;
        counter = counter + 1;
    end;
    Displacementmatrix2 = Localstiffnessmatrix2\Forcematrix;
    Forcematrix3 = Localstiffnessmatrix*Displacementmatrix2;
    Forcematrix3 = Forcematrix3 - Forcematrix2; 
end;