%---------------------------------------------
% Plot initial shape of the beam
%
ux0=XORD';
vy0=zeros(1,NumNP);
plot(ux0,vy0,'-o')
hold on

% -----------------------------------------------------
% plot deformation of the beam
% 
N=50; % N interpolation points in every beam element.

def=[]; % Initial matrix 
i=1; % the 1st beam element
for n=1:N+1 % N+1 interpolation point between the 1st TWO node-points
    x=L/N*(n-1);% local coordinate
    L=abs(XORD(NP(i,2))-XORD(NP(i,1))); % element length
    % get the deflection for the interpolation points in the element
    N1_x = (2*x^3-3*x^2*L+L^3)/L^3; % Shape function
    N2_x = (L*x^3-2*L^2*x^2+x*L^3)/L^3;
    N3_x = (-2*x^3+3*L*x^2)/L^3;
    N4_x = (L*x^3-L^2*x^2)/L^3;
    u1=U(2*NP(i,1)-1);phe1=U(2*NP(i,1));u2=U(2*NP(i,2)-1);phe2=U(2*NP(i,2));
    v_def = [N1_x N2_x N3_x N4_x]*[u1;phe1;u2;phe2];
    def=[def;x v_def];
end
ux=def(:,1)'; % x-coordinate
vy=def(:,2)'; % deflection

for i=2:NumEL % other elements
    for n=1:N % N interpolation point
        x=L/N*n; % local coordinate
        L=abs(XORD(NP(i,2))-XORD(NP(i,1)));% element length
        % get the deflection for the interpolation points in the element
        N1_x = (2*x^3-3*x^2*L+L^3)/L^3; % Shape function
        N2_x = (L*x^3-2*L^2*x^2+x*L^3)/L^3;
        N3_x = (-2*x^3+3*L*x^2)/L^3;
        N4_x = (L*x^3-L^2*x^2)/L^3;
        u1=U(2*NP(i,1)-1);phe1=U(2*NP(i,1));u2=U(2*NP(i,2)-1);phe2=U(2*NP(i,2));
        v_def = [N1_x N2_x N3_x N4_x]*[u1;phe1;u2;phe2];
        def=[def;x+XORD(NP(i,1)) v_def];
    end
    ux=[def(:,1)'];
    vy=[def(:,2)'];
end

% figure,
plot(ux,vy,'r'); grid on
hold on