clear, clc

%INPUT

Lx = 4;

Ly = 2;

E = 200e6;

v = 0.33;

t = 0.01;

q = -5000;

a = Lx/2;

b = Ly/2;

Em = E/(1-v^2)*[1 v 0; v 1 0; 0 0 (1-v)/2]

%SHAPE FUNCTION

syms x y real

N1 = 1-x/a-y/b+x*y/a/b

N2 = x/a-x*y/a/b

N3 = x*y/a/b

N4 = y/b-x*y/a/b

N = [N1 0 N2 0 N3 0 N4 0; 0 N1 0 N2 0 N3 0 N4]

dN1dx = diff(N1,x)

dN1dy = diff(N1,y)

dN2dx = diff(N2,x)

dN2dy = diff(N2,y)

dN3dx = diff(N3,x)

dN3dy = diff(N3,y)

dN4dx = diff(N4,x)

dN4dy = diff(N4,y)

B1 = [dN1dx 0; 0 dN1dy; dN1dy dN1dx]

B2 = [dN2dx 0; 0 dN2dy; dN2dy dN2dx]

B3 = [dN3dx 0; 0 dN3dy; dN3dy dN3dx]

B4 = [dN4dx 0; 0 dN4dy; dN4dy dN4dx]

B = [B1 B2 B3 B4]

%MATRIX k

k = int(int(B’*Em*B,x,0,a),y,0,b)

k = eval(k)

r1 = [0 0 0 0 0 0 0 0]’

r2 = [0 0 0 q*b/2 0 q*b/2 0 0]’

r3 = [0 0 0 0 0 0 0 0]’

r4 = [0 0 0 q*b/2 0 q*b/2 0 0]’

%ASSEMBLE K

K = zeros(18,18)

R = zeros(18,1)

index = [1 2 3 4 9 10 7 8]% element 1

K(index,index) = K(index,index) + k

R(index) = R(index) + r1

index = [3 4 5 6 11 12 9 10]% element 2

K(index,index) = K(index,index) + k

R(index) = R(index) + r2

index = [7 8 9 10 15 16 13 14]% element 3

K(index,index) = K(index,index) + k

R(index) = R(index) + r3

index = [9 10 11 12 17 18 15 16]% element 4

K(index,index) = K(index,index) + k

R(index) = R(index) + r4

%IMPOSE BC

BC = [1 2 7 8 13 14]

K1 = K

K1(BC,:) = []

K1(:,BC) = []

R1 = R

R1(BC) = []

%SOLVE

D1 = t*K1\R1

maxD = max(abs(D1))

D = zeros(18,1)

NBC = 1:18; NBC(BC) = [];

D(NBC) = D1