close all; clear all; clc

qm = 800000;
xL = 1.2;
a = xL/2;
b = 0.1;
h0 = 0.3;
E = 69*10^9;
rho = 0.33;
L = 1.2;


elements = 2;
nodes = 4*elements;
h = L/elements;
xL= L/elements;
x=0:h:L;


xv=0:0.01:1.2;

N1 = 1-3.*(x.^2./xL.^2) + 2.*(x./xL).^3;
N2 = x.*(1-2.*x./xL + (x./xL).^2);
N3 = (x./xL).^2.*(3-2.*x./xL);
N4 = x.^2./xL.*(x./xL-1);

N11 = 1-3.*(xv.^2./L.^2) + 2.*(xv./L).^3;
N22 = xv.*(1-2.*xv./L + (xv./L).^2);
N33 = (xv./L).^2.*(3-2.*xv./L);
N44 = xv.^2./L.*(xv./L-1);
figure()
plot(xv,N11); hold on; plot(xv,N22); hold on; plot(xv,N33); hold on; plot(xv,N44);

dN1 = -6.*xv./xL^2 + 6.*((xv.^2)./xL^3);
dN2 = 1- 4.*xv./xL + 3.*xv.^2./xL^2;
dN3 = 6.*xv./xL^2 - 6.*xv.^2./xL;
dN4 = 3.*xv.^2./xL^2 - 2.*xv./xL;

d2N1 = ((-6./xL^2) + (12.*xv./(xL^3)));
d2N2 = ((-4/xL) + (6.*xv./(xL^2)));
d2N3 = ((6/xL^2) - (12.*xv./xL^3));
d2N4 = ((6.*xv./xL^2) - (2/xL));
figure()
plot(xv,d2N1); hold on; plot(xv,d2N2); hold on; plot(xv,d2N3); hold on; plot(xv,d2N4);

figure()
plot(xv,dN1); hold on; plot(xv,dN2); hold on; plot(xv,dN3); hold on; plot(xv,dN4);

K = zeros(4 + (elements-1)*3); Ku = zeros(4 + (elements-1)*3,1);
i=1; n=i;
while n<length(x)

    Ku1 = @(x) (1-3.*(x.^2./xL.^2) + 2.*(x./xL).^3).*qm.*sin(2*pi*3.*x./(2*xL));
    Ku2 = @(x) (x.*(1-2.*x./xL + (x./xL).^2)).*qm.*sin(2*pi*3.*x./(2*xL));
    Ku3 = @(x) ((x./xL).^2.*(3-2.*x./xL)).*qm.*sin(2*pi*3.*x./(2*xL));
    Ku4 = @(x) (x.^2./xL.*(x./xL-1)).*qm.*sin(2*pi*3.*x./(2*xL));

    K11 = @(x) ((-6./xL^2) + (12.*x./(xL^3))).*E.*((b.*(h0.*exp(-x./xL)).^3)./12).*((-6./xL^2) + (12.*x./(xL^3)));
    K12 = @(x) ((-6./xL^2) + (12.*x./(xL^3))).*E.*(b.*(h0.*exp(-x./xL)).^3)./12.*((-4/xL) + (6.*x./(xL^2)));
    K13 = @(x) ((-6./xL^2) + (12.*x./(xL^3))).*E.*(b.*(h0.*exp(-x./xL)).^3)./12.*((6/xL^2) - (12.*x./xL^3));
    K14 = @(x) ((-6./xL^2) + (12.*x./(xL^3))).*E.*(b.*(h0.*exp(-x./xL)).^3)./12.*((6.*x./xL^2) - (2/xL));
    K22 = @(x) ((-4/xL) + (6.*x./(xL^2))).*E.*(b.*(h0.*exp(-x./xL)).^3)./12.*((-4/xL) + (6.*x./(xL^2)));
    K23 = @(x) ((-4/xL) + (6.*x./(xL^2))).*E.*(b.*(h0.*exp(-x./xL)).^3)./12.*((6/xL^2) - (12.*x./xL^3));
    K24 = @(x) ((-4/xL) + (6.*x./(xL^2))).*E.*(b.*(h0.*exp(-x./xL)).^3)./12.*((6.*x./xL^2) - (2/xL));
    K33 = @(x) ((6/xL^2) - (12.*x./xL^3)).*E.*(b.*(h0.*exp(-x./xL)).^3)./12.*((6/xL^2) - (12.*x./xL^3));
    K34 = @(x) ((6/xL^2) - (12.*x./xL^3)).*E.*(b.*(h0.*exp(-x./xL)).^3)./12.*((6.*x./xL^2) - (2/xL));
    K44 = @(x) ((6.*x./xL^2) - (2/xL)).*E.*(b.*(h0.*exp(-x./xL)).^3)./12.*((6.*x./xL^2) - (2/xL));

    X1 = x(n);
    X2 = x(n+1);

    K(i,i) = K(i+i) + integral(K11, X1, X2);
    K(i,i+1) = integral(K12, X1, X2); K(i+1,i) = K(i,i+1);
    K(i,i+2) = integral(K13, X1, X2); K(i+2,i) = K(i,i+2);
    K(i,i+3) = integral(K14, X1, X2); K(i+3,i) = K(i,i+3);
    K(i+1,i+1) = integral(K22, X1, X2);
    K(i+1,i+2) = integral(K23, X1, X2); K(i+2,i+1) = K(i+1,i+2);
    K(i+1,i+3) = integral(K24, X1, X2); K(i+3,i+1) = K(i+1,i+3);
    K(i+2,i+2) = integral(K33, X1, X2);
    K(i+2,i+3) = integral(K34, X1, X2); K(i+3,i+2) = K(i+2,i+3);
    K(i+3,i+3) = integral(K44, X1, X2);

    if x(n)>=2*L/3
        Ku(i,1) = K(i+1) + integral(Ku1, X1, X2);
        Ku(i+1,1) = integral(Ku2, X1, X2);
        Ku(i+2,1) = integral(Ku3, X1, X2);
        Ku(i+3,1) = integral(Ku4, X1, X2);
    end

    i= i+3; n= n+1;
end
K(1,1) = K(1,1)*10^10; K(2,2) = K(1,1);

    w=K\Ku;
    xx = 0:L/length(w):L-L/length(w);
    figure()
    plot(xx,w)