clear all; close all; clc
%
%syms  N1(r,s) N2(r,s) N3(r,s) N4(r,s)
syms r s
N1=1/4*(1+r)*(1+s);
N2=1/4*(1-r)*(1+s);
N3=1/4*(1-r)*(1-s); 
N4=1/4*(1+r)*(1-s);
Jacob=[20,0;0,20];
Poisson = 0.3;
Elast=1;
Force=[0; 0; 0; 0; -1; 0; 0; 0];
D_Matrix=Elast/(1-Poisson^2)*[1, Poisson, 0; Poisson, 1, 0; 0, 0, (1-Poisson)/2];
dN1dx=20/det(Jacob)*diff(N1,r);
dN1dy=20/det(Jacob)*diff(N1,s);
dN2dx=20/det(Jacob)*diff(N2,r);
dN2dy=20/det(Jacob)*diff(N2,s);
dN3dx=20/det(Jacob)*diff(N3,r);
dN3dy=20/det(Jacob)*diff(N3,s);
dN4dx=20/det(Jacob)*diff(N4,r);
dN4dy=20/det(Jacob)*diff(N4,s);

B_Matrix = [dN1dx, dN2dx, dN3dx, dN4dx, 0, 0, 0, 0; ...
            0, 0, 0, 0, dN1dy, dN2dy, dN3dy, dN4dy; ...
            dN1dy, dN2dy, dN3dy, dN4dy, dN1dx, dN2dx, dN3dx, dN4dx];
fun = B_Matrix'*D_Matrix*B_Matrix*det(Jacob);
K = int(int(fun,s,-1,1),r,-1,1);
K_comp=K;


K(2,:)=[];
K(:,2)=[];
Force(2)=[];
K(2,:)=[];
K(:,2)=[];
Force(2)=[];
K(5,:)=[];
K(:,5)=[];
Force(5)=[];
K

Force
Displacement = double(inv(K)*Force)
U=[Displacement(1); 0; 0; Displacement(2); Displacement(3); Displacement(4); 0;  Displacement(5)]
Reactions = double(K_comp*U)