clear() // clear all variables
rho_m = 3330;
rho_fill = 2400;
g = 9.81;
Pb = -1.5e12;
v = 0.25;
E = 1e11;

Te = [5000;10000;25000]; // Elastic thickness (m)
D = E*Te^3/(12*(1-v^2)); // flexural rigidity
lamda = ((rho_m-rho_fill)*g./(4*D)).^0.25; // lamda = [(rho_m - rho_infill)*g/4*D]^0.25
lamda_k(:,1) = lamda(:,1)*1000;
dist = [0:5:200]; // x-limits
flex(:,1) = dist'; // defines the x-axis format and limits
flex(:,2) = 2*Pb*lamda(1,1)/((rho_m-rho_fill)*g)*exp(-lamda_k(1,1)*flex(:,1)).*cos(lamda_k(1,1)*flex(:,1));
flex(:,3) = 2*Pb*lamda(2,1)/((rho_m-rho_fill)*g)*exp(-lamda_k(2,1)*flex(:,1)).*cos(lamda_k(2,1)*flex(:,1));
flex(:,4) = 2*Pb*lamda(3,1)/((rho_m-rho_fill)*g)*exp(-lamda_k(3,1)*flex(:,1)).*cos(lamda_k(3,1)*flex(:,1));
scf(1) // define Graphic window 1
clf(1)
set(gca(),"auto_clear","off")
set(gca(),"grid",[0,0])
plot(flex(:,1),flex(:,2),flex(:,1),flex(:,3),flex(:,1),flex(:,4))
legend('$T_e\ = 0\ km$','$T_e\ = 5\ km$', '$T_e\ = 10\ km$', '$T_e\ = 25\ km$',4);
xlabel('$Distance\ (km)$','fontsize',5.5)
ylabel('$Flexure\ (m)$','fontsize',5.5)
title('$Flexure\ for\ a\ Broken\ Plate\ Model$','fontsize',5.5)