diff --git a/README.md b/README.md
index bf372d9..76b1062 100644
--- a/README.md
+++ b/README.md
@@ -102,3 +102,73 @@ 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 = pw-se;
 ```
+
+# Part F
+```matlab
+ n=[3,20:5:40];
+ P=0.001; %MPa
+ T = zeros(1,length(n));
+ ea = zeros(1,length(n));
+ for i = 1:length(n)
+     [T(i), ea(i)] = tension_sol(P,n(i));
+ end
+ ```
+ ```matlab
+ function [T,ea] = tension_sol(P,n)
+y =@(T) SE_diff(T,P,n);
+[T,fx,ea,iter]=bisect(y,.01,1);
+```
+```matlab
+function [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,varargin)
+% bisect: root location zeroes
+% [root,fx,ea,iter]=bisect(func,xl,xu,es,maxit,p1,p2,...):
+% uses bisection method to find the root of func
+% input:
+% func = name of function
+% xl, xu = lower and upper guesses
+% es = desired relative error (default = 0.0001%)
+% maxit = maximum allowable iterations (default = 50)
+% p1,p2,... = additional parameters used by func
+% output:
+% root = real root
+% fx = function value at root
+% ea = approximate relative error (%)
+% iter = number of iterations
+if nargin<3,error('at least 3 input arguments required'),end
+test = func(xl,varargin{:})*func(xu,varargin{:});
+if test>0,error('no sign change'),end
+if nargin<4||isempty(es), es=0.0001;end
+if nargin<5||isempty(maxit), maxit=50;end
+iter = 0; xr = xl; ea = 100;
+while (1)
+  xrold = xr;
+  xr = (xl + xu)/2;
+  iter = iter + 1;
+  if xr ~= 0,ea = abs((xr - xrold)/xr) * 100;end
+  test = func(xl,varargin{:})*func(xr,varargin{:});
+  if test < 0
+    xu = xr;
+  elseif test > 0
+    xl = xr;
+  else
+    ea = 0;
+  end
+  if ea <= es || iter >= maxit,break,end
+end
+root = xr; fx = func(xr, varargin{:});
+```
+```matlab
+function re = Rel_error (T)
+re = zeros(1,length(T)-1);
+for i = 2:length(T)
+    re(i-1)= abs(T(i)-T(i-1))/T(i-1);
+end
+```
+|number of nodes|Tension (uN/um)| rel. error|
+|---|---|---|
+|3 |0.0489 |n/a|
+|20|0.0599|22.6%|
+|25|0.0601|0.27%|
+|30|0.0602|0.15%|
+|35|0.0602|0.09%|
+|40|0.0603|0.06%|