diff --git a/drawFunc.py b/drawFunc.py
index 02e9777..1b1a67b 100644
--- a/drawFunc.py
+++ b/drawFunc.py
@@ -16,5 +16,6 @@ def draw(f_x):
 
 # take this pencil
     levels = [0.0, 1, 2, 5, 20, 30, 50, 80, 100, 200, 400]
-    cs = plt.contour(X, Y, Z, levels, colors='k')
-    plt.clabel(cs, fmt='%.1f', inline=1)
+    #cs = plt.contour(X, Y, Z, levels, colors='k')
+    cs = plt.contourf(X, Y, Z, 20, cmap=plt.cm.YlGnBu)
+    #plt.clabel(cs, fmt='%.1f', inline=1)
diff --git a/func.py b/func.py
index cad1fdb..6a32d9e 100644
--- a/func.py
+++ b/func.py
@@ -31,7 +31,7 @@ def __init__(self):
         pass
 
     def f_x(self, x):
-        f = 10 * (x[1] - (x[0])**2)**2 + (1-x[0])**2
+        f = 10.0 * (x[1] - (x[0])**2)**2 + (1-x[0])**2
         return f[0, 0]
 
     def g_x(self, x):
@@ -40,7 +40,7 @@ def g_x(self, x):
         return np.matrix([[g0[0,0]], [g1[0,0]]])
 
     def G_x(self, x):
-        G00 = 120.0 * x[0]**2 - 40.0*x[1] + 2
+        G00 = 120.0 * x[0]**2 - 40.0 * x[1] + 2.0
         G01 = -40.0*x[0]
         G10 = -40.0*x[0]
         G11 = 20.0
diff --git a/main.py b/main.py
index ace3dda..106a12f 100644
--- a/main.py
+++ b/main.py
@@ -34,17 +34,17 @@ def rosen():
     plt.figure()
 
     drawFunc.draw(obj.f_x)
-    start_point = np.matrix([[-1.5], [0.0]])
-    plt.annotate('Start', xy=(-1.5, 0.0), xytext=(-1.8, 0.5),
+    start_point = np.matrix([[-1.5], [-4.0]])
+    plt.annotate('Start', xy=(-1.5, -4.0), xytext=(-1.8, -3.5),
             arrowprops=dict(facecolor='black', shrink=0.02, frac=0.5))
     plt.annotate('Optimal', xy=(1, 1), xytext=(1.1, 1.4),
             arrowprops=dict(facecolor='black', shrink=0.02, frac=0.5))
 
-    v = 0.0
-    str1 = "Newton"
-    track = opt.newton(start_point, obj, v)
-    p1, = plt.plot(track[0,:].tolist()[0], track[1,:].tolist()[0],
-            'o-', linewidth=2.0)
+    #v = 0.0
+    #str1 = "Newton"
+    #track = opt.newton(start_point, obj, v)
+    #p1, = plt.plot(track[0,:].tolist()[0], track[1,:].tolist()[0],
+    #        'o-', linewidth=2.0)
 
     v = 0.1
     str2 = "Modified_Newton %.2f"%v
@@ -68,18 +68,15 @@ def rosen():
     p5, = plt.plot(track[0,:].tolist()[0], track[1,:].tolist()[0],
             'o-', linewidth=2.0)
 
-    str6 = "BFGS_trivial"
-    track = opt.quasi_newton(start_point, obj, 150, 1)
+    str6 = "BFGS_Armijo"
+    track = opt.BFGS(start_point, obj, 350, 0)
     p6, = plt.plot(track[0,:].tolist()[0], track[1,:].tolist()[0],
             'o-', linewidth=2.0)
 
-    str7 = "BFGS_trivial"
-    track = opt.quasi_newton(start_point, obj, 150, 0)
-    p7, = plt.plot(track[0,:].tolist()[0], track[1,:].tolist()[0],
-            'o-', linewidth=2.0)
-
-    plt.legend([p1, p2, p3, p4, p5, p6, p7],\
-            [str1, str2, str3, str4, str5, str6, str7])
+    #plt.legend([p1, p2, p3, p4, p5, p6, p7],\
+    #        [str1, str2, str3, str4, str5, str6, str7])
+    plt.legend([p2, p3, p4, p5, p6],\
+            [str2, str3, str4, str5, str6])
     plt.show()
 
 #quadratic()
diff --git a/optMethods.py b/optMethods.py
index ec31bdb..0950495 100644
--- a/optMethods.py
+++ b/optMethods.py
@@ -1,4 +1,5 @@
-import numpy as np
+from numpy import *
+import pdb
 
 def newton(start_point, obj_fun, modified=0, iterations=5):
 # with second order information
@@ -9,15 +10,15 @@ def newton(start_point, obj_fun, modified=0, iterations=5):
     while k < iterations:
 
         if modified > 0:
-            vI = modified * np.matrix([[1, 0], [0, 1]])
-            G_ = np.linalg.inv(obj_fun.G_x(x) + vI) # inverse?
+            vI = modified * matrix([[1, 0], [0, 1]])
+            G_ = linalg.inv(obj_fun.G_x(x) + vI) # inverse?
         else:
-            G_ = np.linalg.inv(obj_fun.G_x(x)) # inverse?
+            G_ = linalg.inv(obj_fun.G_x(x)) # inverse?
         g_ = obj_fun.g_x(x)
-        delta = -1.0 * np.dot(G_, g_)
+        delta = -1.0 * dot(G_, g_)
         x = x + delta
 
-        track = np.concatenate((track, x), axis=1)
+        track = concatenate((track, x), axis=1)
         k += 1
 
     return track
@@ -28,48 +29,47 @@ def BFGS(start_point, obj_fun, iteration=10, interpolation=0):
     track = x
 
     k = 0
-    H = np.matrix([[.1, 0], [0, .1]])
+    #H = dot(obj_fun.g_x(x), obj_fun.g_x(x).T)
+    H = matrix([[.10, .0], [.0, .10]])
+    gamma = 1.0
 
-    while k < iteration:
+    while k < iteration and linalg.norm(gamma) > 1e-10:
 
-        s = -1.0 * np.dot(H, obj_fun.g_x(x))
+        p = -dot(H, obj_fun.g_x(x))
 
         if interpolation > 0:
-            alpha_k = _backtracking_line_search(obj_fun, x, s)
+            alpha_k = _backtracking_line_search(obj_fun, x, p)
         else:
-            alpha_k = _armijo_line_search(obj_fun, x, s)
+            alpha_k = _armijo_line_search(obj_fun, x, p)
 
-        delta = alpha_k * s
-        x_k_1 = x + delta
+        s = alpha_k * p
+        g_k = obj_fun.g_x(x)
+        x = x + alpha_k * p
+        g_k_1 = obj_fun.g_x(x)
 
-        gamma = obj_fun.g_x(x_k_1) - obj_fun.g_x(x)
-        delta_T_gamma = np.dot(delta.T, gamma)
-        gamma_T_H_gamma = np.dot(np.dot(gamma.T, H), gamma)
-        delta_delta_T = np.dot(delta, delta.T)
-
-        H = H + (1 + gamma_T_H_gamma/delta_T_gamma) *\
-                (delta_delta_T / delta_T_gamma)
+        y = g_k_1 - g_k
 
-        delta_gamma_H = np.dot(np.dot(delta,delta.T), H)
-        H_gamma_delta_T = np.dot(np.dot(H, gamma), delta.T)
-        H = H - (delta_gamma_H + H_gamma_delta_T) / delta_T_gamma
+        z = dot(H, y)
+        sTy = dot(s.T, y)
+        H += outer(s, s) * (sTy + dot(y.T, z))[0,0]/(sTy**2) \
+                - (outer(z, s) + outer(s, z))/sTy
 
-        track = np.concatenate((track, x), axis=1)
-        x = x_k_1
+        track = concatenate((track, x), axis=1)
         k += 1
 
     return track
+
 def quasi_newton(start_point, obj_fun, iteration=10, interpolation=0):
 # with first order information
     x = start_point
     track = x
 
     k = 0
-    H = np.matrix([[.1, 0], [0, .1]])
+    H = matrix([[.01, 0], [0, .01]])
 
     while k < iteration:
 
-        s = -1.0 * np.dot(H, obj_fun.g_x(x))
+        s = -1.0 * dot(H, obj_fun.g_x(x))
 
         if interpolation > 0:
             alpha_k = _backtracking_line_search(obj_fun, x, s)
@@ -80,14 +80,14 @@ def quasi_newton(start_point, obj_fun, iteration=10, interpolation=0):
         x_k_1 = x + delta
 
         gamma = obj_fun.g_x(x_k_1) - obj_fun.g_x(x)
-        u = delta - np.dot(H, gamma)
+        u = delta - dot(H, gamma)
 
-        scale_a = np.dot(u.T, gamma)
+        scale_a = dot(u.T, gamma)
         if scale_a == 0: # :(
             scale_a = 0.000001
-        H = H + np.dot(u, u.T) / scale_a
+        H = H + dot(u, u.T) / scale_a
 
-        track = np.concatenate((track, x), axis=1)
+        track = concatenate((track, x), axis=1)
         x = x_k_1
         k += 1
 
@@ -111,9 +111,8 @@ def _backtracking_line_search(obj_fun, x, s, c=1e-4):
     return alpha
 
 
-def _armijo_line_search(obj_fun, x, s, c=1e-4):
+def _armijo_line_search(obj_fun, x, s, c=1e-4, alpha0=1.0):
 # adapted from scipy.optimisation
-    alpha0 = 1.0
     amin = 0.0
     f_alpha = obj_fun.about_alpha(x, s)
     g_alpha = obj_fun.about_alpha_prime(x, s)
@@ -137,7 +136,7 @@ def _armijo_line_search(obj_fun, x, s, c=1e-4):
             alpha1**3 * (f_alpha(alpha0) - f_alpha(0) - g_alpha(0) * alpha0)
         b = b / factor
 
-        alpha2 = (-b + np.sqrt(abs(b**2 - 3 * a * g_alpha(0)))) / (3.0 * a)
+        alpha2 = (-b + sqrt(abs(b**2 - 3 * a * g_alpha(0)))) / (3.0 * a)
         if(f_alpha(alpha2) <= f_alpha(0) + c * alpha2 * g_alpha(0)):
             return alpha2