diff --git a/PS1/q1b.m b/PS1/q1b.m
index 496bb32..1c06119 100644
--- a/PS1/q1b.m
+++ b/PS1/q1b.m
@@ -3,18 +3,18 @@ clear all;
 %% Constants
 
 % Physical constants
-hbar = 1.052e-34;
-q = 1.602e-19;
-epsilon_0 = 8.854e-12;
-epsilon_r = 4;
-mstar = 0.25 * 9.11e-31;
+%hbar = 1.052e-34;
+%q = 1.602e-19;
+%epsilon_0 = 8.854e-12;
+%epsilon_r = 4;
+%mstar = 0.25 * 9.11e-31;
 
 % Single-charge coupling energy (eV)
 U_0 = 0.25;
 % (eV)
 kBT = 0.025;
 % Contact coupling coefficients (eV)
-gamma_1 = 0.0005;
+gamma_1 = 0.005;
 gamma_2 = gamma_1;
 gamma_sum = gamma_1 + gamma_2;
 % Capacitive gate coefficient
@@ -35,13 +35,11 @@ cal_E = 0.2;
 
 % Lorentzian density of states, normalized so the integral is 1
 D = (gamma_sum / (2*pi)) ./ ( (E-cal_E).^2 + (gamma_sum/2).^2);
-
+D = D ./ (dE*sum(D));
 
 % Reference no. of electrons in channel
 N_0 = 0;
 
-fermi(-0.25, -0.2, kb_T)
-
 voltages = linspace(0, 1, 101);
 
 % Terminal Voltages
@@ -52,24 +50,28 @@ for n = 1:length(voltages)
     % Set varying drain voltage
     V_D = voltages(n);
 
+    % Shifted energy levels of the contacts
     mu_1 = mu - V_S;
     mu_2 = mu - V_D;
 
     % Laplace potential, does not change as solution is found (eV)
-    U_L = -q * ((C_S*V_S + C_G*V_G + C_D*V_D) / C_E);
+    U_L = - (a_G*V_G - a_D*V_D - a_S*V_S);
 
     % Poisson potential must change, assume 0 initially (eV)
     U_P = 0;
 
+    % Assume large rate of change
     dU_P = 1;
+
+    % Run until we get close enough to the answer
     while dU_P > 1e-6
         % source Fermi function
-        f_1 = 1 / (1 + exp((E + U_L + U_P - mu_1) / kBT));
+        f_1 = 1 ./ (1 + exp((E + U_L + U_P - mu_1) ./ kBT));
         % drain Fermi function
-        f_2 = 1 / (1 + exp((E + U_L + U_P - mu_2) / kBT));
+        f_2 = 1 ./ (1 + exp((E + U_L + U_P - mu_2) ./ kBT));
 
         % Update channel electrons against potential
-        N(n) = dE * sum( ((gamma_1/gamma_sum) .* f1 + (gamma_2/gamma_sum) .* f2) .* D);
+        N(n) = dE * sum( ((gamma_1/gamma_sum) .* f_1 + (gamma_2/gamma_sum) .* f_2) .* D);
 
         % Re-update Poisson portion of potential
         tmpU_P = U_0 * ( N(n) - N_0);
@@ -80,7 +82,9 @@ for n = 1:length(voltages)
         % U_P = U_P + 0.1 * (tmpU_P - U_P)
     end
 
-    I(n) = q * (q/hbar) * (gamma_1 * gamma_1 / gamma_sum) * dE * sum((f1-f2).*D);
+    % Calculate current based on solved potential.
+    % Note: f1 is dependent on changes in U but has been updated prior in the loop
+    I(n) = q * (q/hbar) * (gamma_1 * gamma_1 / gamma_sum) * dE * sum((f_1-f_2).*D);
 
 end