CODE
clc
clear all
format long

EI=2400;
L=1;
qo=-600;

% elements
n=[1,2,4,10,30,50,100];

for division=7

    K=zeros(2*(n(division)+1));
    X=Linspace(0,1,n(division));
    K_element=[12*EI/(L/n(division))^3,6*EI/(L/n(division))^2, -12*EI/(1/n(division))^3,6*EI/(1/n(division))^2;6*EI/(1/n(division))^2,4*EI/(L/n(division)),-6*EI/(1/n(division))^2,2*EI/(L/n(division));-12*EI/(L/n(division))^3,-6*EI/(1/n(division))^2,12*EI/(L/n(division))^3,-6*EI/(1/n(division))^2;6*EI/(1/n(division))^2,2*EI/(L/n(division)),-6*EI/(1/n(division))^2,4*EI/(L/n(division))];


    % Boundary values of K
    for i=1:2

        for j=1:4

            K(i,j)=K_element(i,j);
            K(i + 2*n(division),j + 2*n(division)-2)=K_element(i+2,j);
        end
    end

    % (2i-1)th values for K displacement
    for i=3:2*n(division)
        for j=1:2*n(division)+2
            if mod(i,2)~=0
                if i==j+2
                    K(i,j)=K_element(3,1);
                else if i==j+1
                        K(i,j)=K_element(3,2);
                    else if i==j
                            K(i,j)=K_element(1,1)+K_element(3,3);
                        else if i==j-1
                                K(i,j)=K_element(3,4)+K_element(1,2);
                            else if i==j-2
                                    K(i,j)=K_element(1,3);
                                else if i==j-3
                                        K(i,j)=K_element(1,4);
                                    else
                                        K(i,j)=0;
                                    end
                                end
                            end
                        end
                    end
                end
            end

            % 2ith values for angle
            if mod(i,2)==0
                if i==j+3
                    K(i,j)=K_element(4,1);
               else if i==j+2
                       K(i,j)=K_element(4,2);
                   else if i==j+1
                           K(i,j)=K_element(4,3)+ K_element(2,1);
                       else if i==j
                               K(i,j)=K_element(4,4)+K_element(2,2);
                           else if i==j-1
                                   K(i,j)=K_element(2,3);
                               else if i==j-2
                                       K(i,j)=K_element(2,4);
                                   else
                                       K(i,j)=0;
                                   end
                               end
                           end
                       end
                   end
               end
           end
        end
    end

    K_backup=K;
    % Boundary values for force and moment
    F(1)=qo*L/(2*n(division));
    F(2)=qo*(L/n(division))^2/12;
    F(2*n(division)+1)=qo*L/(2*n(division));
    F(2*n(division)+2)=-1*qo*(L/n(division))^2/12;

    % (2i-1)th and 2ith values for force and moment
    for i=2:n(division)
        F(2*i-1)=qo*L/n(division);
        F(2*i)=0;
    end





    %K(:,2*n(division)+2)=[];
    K(:,2*n(division)+1)=[];
    %K(2*n(division)+2,:)=[];
    K(2*n(division)+1,:)=[];

    %K(:,n(division)+1)=[];
    %K(n(division)+1,:)=[];

    %K(:,2)=[];
    K(:,1)=[];
    %K(2,:)=[];
    K(1,:)=[];


    %F(2*n(division)+2)=[];
    F(2*n(division)+1)=[];
    %F(n(division)+1)=[];
    %F(2)=[];
    F(1)=[];

    d=K\F';

    % simple support
    d1(1)=0;

    % simple support
    %d1(n(division))=0;

    for i=1:n(division)-1
        d1(i)=d(2*i);
    end

    d1(n(division))=0;

    % simple support
    % d1(n(division)/2+1)=0;

    % for i=(n(division)/2):n(division)-2
    %    d1(i+2)=d(2*i+1);
    % end

    % simple support
    % d1(n(division)+1)=0;

    plot(X,d1)

end