diff --git a/HW3_Script.m b/HW3_Script.m
index a1b4018..b277545 100644
--- a/HW3_Script.m
+++ b/HW3_Script.m
@@ -1,10 +1,10 @@
 %HW_3
 %Script calls functions and outputs tables and plots relevant to HW3.
 
- [root(1),fx(1),ea(1),iter(1)]=falsepos(@(x) projectile(x),1,50,.00001);
- [root(2),ea(2),iter(2)]=mod_secant(@(x) projectile(x),1,50,.00001);
- [root(3), ea(3), iter(3)]= newtraph(@(x) projectile(x),@(x) dprojectile_dtheta(x),50,.00001,1000);
- [root(4),fx(4),ea(4),iter(4)]=bisect(@(x) projectile(x),0,50,.00001);
+ [root1,fx1,ea1,iter1]=falsepos(@(x) projectile(x),1,50,.00001);
+ [root2,ea2,iter2]=mod_secant(@(x) projectile(x),1,50,.00001);
+ [root3, ea3, iter3]= newtraph(@(x) projectile(x),@(x) dprojectile_dtheta(x),50,.00001,1000);
+ [root4,fx4,ea4,iter4]=bisect(@(x) projectile(x),0,50,.00001);
 
 
  %Question 1: part (i)
@@ -19,20 +19,23 @@
 
 %Question 1: part(ii)
 %Plots the graphs that compare methods
-
-    for n = 1:iter(1)
+e_b = 1;
+e_f = 2;
+e_n = 3;
+e_m = 4;
+    for n = 1:iter1
         [r, y, e_b(n), k] = falsepos(@(x) projectile(x),1,50,.00001);
     end
     
-    for n = 1:iter(2)
+    for n = 1:iter2
         [r, y, e_f(n)] = mod_secant(@(x) projectile(x),1,50,.00001);
     end
     
-    for n = 1:iter(3)
+    for n = 1:iter3
         [r, e_n(n), k] = newtraph(@(x) projectile(x),@(x) dprojectile_dtheta(x),50,.00001,1000);
     end
     
-    for n = 1:iter(4)
+    for n = 1:iter4
         [r, e_m(n), k] = bisect(@(x) projectile(x),0,50,.00001);
     end
     
@@ -40,7 +43,7 @@
  %This section sets up the code inorder to plot the data.
     
     setdefaults
-    plot(1:iter(1), e_b, 'g', 1:iter(2), e_f, 'b', 1:iter(3), e_n, 'r', 1:iter(4), e_m, 'y')
+    plot(1:iter1, e_b, 'g', 1:iter2, e_f, 'b', 1:iter3, e_n, 'r', 1:iter4, e_m, 'y')
     title('Approx Err. of Convergent Functions vs Number of Iterations');
     xlabel('Iteration');
     ylabel('Approx Error');
diff --git a/HW4_Script.m b/HW4_Script.m
new file mode 100644
index 0000000..346b95f
--- /dev/null
+++ b/HW4_Script.m
@@ -0,0 +1,28 @@
+%HW_4 Script
+
+setdefaults;
+
+%Part B: using the goldmin function
+theta = 0;
+[x1,fx1,ea1,iter1] = goldmin(@(x) collar_potential_energy(x,theta),-10,10,0.0001);
+
+%Part C: Creating For-loop solving for potential energy
+
+theta = .9*pi/180;  
+thetaVector = linspace(0,pi, 400);
+fx2 = 0;
+
+[x3,fx3,ea3,iter3] = goldmin(@(x) collar_potential_energy(x,theta),-10,10,0.0001);
+
+for n = 1:length(thetaVector)
+    theta = thetaVector(n);
+[x2,fx2(n),ea2,iter2] = goldmin(@(x) collar_potential_energy(x,theta),-10,10,0.0001);
+end
+
+%Part D: Plotting the Data
+
+plot(thetaVector, fx2);
+    title('Plot of Minimum PE for x_C with given Angle');
+    xlabel('Theta (radians)');
+    ylabel('x_C (meters)');
+    
\ No newline at end of file
diff --git a/collar_potential_energy.m b/collar_potential_energy.m
new file mode 100644
index 0000000..37b72dc
--- /dev/null
+++ b/collar_potential_energy.m
@@ -0,0 +1,16 @@
+function [ PE_tot ] = collar_potential_energy( x_C, theta )
+
+% Computes total Potential energy of the collar system given
+% an initial x_C and theta
+
+m = 0.5; % mass in kg
+g = 9.81; % acceleration due to gravity m/s^2
+K = 30;  %spring constant N/m
+
+PE_g = m*x_C*g*sin(theta);
+PE_s = .5*K *((0.5 - sqrt(0.5^2+(0.5-x_C)^2)))^2;
+
+PE_tot = PE_g + PE_s;
+
+end
+
diff --git a/goldmin.m b/goldmin.m
new file mode 100644
index 0000000..fb4ce0b
--- /dev/null
+++ b/goldmin.m
@@ -0,0 +1,36 @@
+function [x,fx,ea,iter]=goldmin(f,xl,xu,es,maxit,varargin)
+% goldmin: minimization golden section search
+%   [x,fx,ea,iter]=goldmin(f,xl,xu,es,maxit,p1,p2,...):
+%     uses golden section search to find the minimum of f
+% input:
+%   f = 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 f
+% output:
+%   x = location of minimum
+%   fx = minimum function value
+%   ea = approximate relative error (%)
+%   iter = number of iterations
+if nargin<3,error('at least 3 input arguments required'),end
+if nargin<4|isempty(es), es=0.0001;end
+if nargin<5|isempty(maxit), maxit=50;end
+phi=(1+sqrt(5))/2;
+iter=0;
+while(1)
+  d = (phi-1)*(xu - xl);
+  x1 = xl + d;
+  x2 = xu - d;
+  if f(x1,varargin{:}) < f(x2,varargin{:})
+    xopt = x1;
+    xl = x2;
+  else
+    xopt = x2;
+    xu = x1;
+  end
+  iter=iter+1;
+  if xopt~=0, ea = (2 - phi) * abs((xu - xl) / xopt) * 100;end
+  if ea <= es | iter >= maxit,break,end
+end
+x=xopt;fx=f(xopt,varargin{:});