Part F

Part F/F_sol.m

@@ -1,4 +1,4 @@
- n=20:5:40;
+ n=[3,20:5:40];
  P=0.001; %MPa
  T = zeros(1,length(n));
  ea = zeros(1,length(n));

Part F/SE_diff.m

@@ -18,9 +18,9 @@ dwdx = zeros(num);
 dwdy = zeros(num);
 for i = 1:num
     for j = 1:num
-        dwdx(i,j) = mean([z(i+1,j)-z(i,j),z(i+1,j+1)-z(i,j+1)]);
-        dwdy(i,j) = mean([z(i,j+1)-z(i,j),z(i+1,j+1)-z(i+1,j)]);
+        dwdx(i,j) = mean([z(i+1,j)-z(i,j),z(i+1,j+1)-z(i,j+1)])./h;
+        dwdy(i,j) = mean([z(i,j+1)-z(i,j),z(i+1,j+1)-z(i+1,j)])./h;
     end
 end
 se = E*t*h^2/(2*(1-v^2))*sum(sum(0.25.*dwdx.^4+.25.*dwdy.^4+0.5.*(dwdx.*dwdy).^2));
-pw_se = abs(pw-se);
+pw_se = pw-se;

Part F/tension_sol.m

@@ -1,3 +1,3 @@
 function [T,ea] = tension_sol(P,n)
 y =@(T) SE_diff(T,P,n);
-[T,fx,ea,iter]=bisect(y,.01,1);
+[T,fx,ea,iter]=bisect(y,.01,1,.1);