diff --git a/fire2d.py b/fire2d.py
index a5a957a1a41067752ebdc5c49610403c15f53130..dafd15d6e24a2798dc9c8ac76b6098b8ff397f58 100644
--- a/fire2d.py
+++ b/fire2d.py
@@ -6,6 +6,7 @@ from numpy.linalg import inv
 import math
 from scipy.integrate import solve_ivp
 from matplotlib import pyplot as plt
+from scipy.optimize import minimize
 
 def sum_to(k):
     return k*(k+1)/2
@@ -102,8 +103,23 @@ def solve(M, y0, T, delta_t):
 
     return sln, i+1
 
+# def squared_change_rate(y, M):
+#     dy_dt = np.matmul(M, y) + np.exp(y)
+#     return np.dot(dy_dt, dy_dt)
+
+# def sln_change(y, M):
+#     sln, i = solve(M, y, 10, 0.1)
+#     change = np.abs(sln[10,:] - y) - np.abs(sln[99,:] - y) - 1e6 * i**2
+#     return (change**2).sum()
+
+# def get_equillibrium(M):
+#     ymin = minimize(lambda y: sln_change(y, M),
+#                     1*np.ones(M.shape[0]), method="Nelder-Mead")
+#     print(ymin)
+#     return ymin.x
+
 if __name__ == "__main__":
-    mesh = Mesh2D(50, 50, 40)
+    mesh = Mesh2D(50, 50, 60)
     # for k in range(mesh.N):
     #     for j in range(mesh.N-k):
     #         print("{:<5}".format(mesh[j,k]), end="")
@@ -112,7 +128,7 @@ if __name__ == "__main__":
     M = construct_matrix(mesh) / Da
 
     ic_low = 0 #3.4
-    ic_high = 10 #3.4
+    ic_high = 0 #3.4
     i = 0
     
     while i < 100:
@@ -125,10 +141,11 @@ if __name__ == "__main__":
         sln, N_t_in_sln = solve(M, initial_conditions, simulation_time, time_step)
 
         print(ic, N_t_in_sln)
+
         if N_t_in_sln < round(simulation_time / time_step):
             ic_high = ic
         else:
-            if i > 50:
+            if i > 0:
                 break
             ic_low = ic
 
@@ -140,7 +157,7 @@ if __name__ == "__main__":
     P = mesh.get_plot_matrix(y_final)
 
     temp_min = 0.0
-    temp_max = 6
+    temp_max = 0.5
 
      # Determine subplot grid
     ncols = 4