diff --git a/src/NewtonKrylov.jl b/src/NewtonKrylov.jl
index ceadb42..2bb7401 100644
--- a/src/NewtonKrylov.jl
+++ b/src/NewtonKrylov.jl
@@ -121,6 +121,87 @@ function (F::EisenstatWalker)(η, tol, n_res, n_res_prior)
 end
 inital(F::EisenstatWalker) = F.η_max
 
+function update_lambda(iarm, armfix, lambda, lamc, ff0, ffc, ffm)
+    if iarm == 0 || armfix == true
+        lambda = lambda * 0.5
+    else
+        lamm = lamc
+        lamc = lambda
+        lambda = parab3p(lamc, lamm, ff0, ffc, ffm)
+    end
+    return lambda
+end
+
+"""
+parab3p(lambdac, lambdam, ff0, ffc, ffm)
+
+Three point parabolic line search.
+
+input:\n
+       lambdac = current steplength
+       lambdam = previous steplength
+       ff0 = value of || F(x_c) ||^2
+       ffc = value of || F(x_c + lambdac d) ||^2
+       ffm = value of || F(x_c + lambdam d) ||^2
+
+output:\n
+       lambdap = new value of lambda
+
+internal parameters:\n
+       sigma0 = .1, sigma1=.5, safeguarding bounds for the linesearch
+
+You get here if cutting the steplength in half doesn't get you
+sufficient decrease. Now you have three points and can build a parabolic
+model. I do not like cubic models because they either need four points
+or a derivative. 
+
+So let's think about how this works. I cheat a bit and check the model
+for negative curvature, which I don't want to see.
+
+ The polynomial is
+
+ p(lambda) = ff0 + (c1 lambda + c2 lambda^2)/d1
+
+ d1 = (lambdac - lambdam)*lambdac*lambdam < 0
+ So if c2 > 0 we have negative curvature and default to
+      lambdap = sigma0 * lambda
+ The logic is that negative curvature is telling us that
+ the polynomial model is not helping much, so it looks better
+ to take the smallest possible step. This is not what I did in the
+ matlab code because I did it wrong. I have sinced fixed it.
+
+ So (Students, listen up!) if c2 < 0 then all we gotta do is minimize
+ (c1 lambda + c2 lambda^2)/d1 over [.1* lambdac, .5*lambdac]
+ This means to MAXIMIZE c1 lambda + c2 lambda^2 becase d1 < 0.
+ So I find the zero of the derivative and check the endpoints.
+
+"""
+function parab3p(lambdac, lambdam, ff0, ffc, ffm)
+    #
+    # internal parameters
+    #
+    sigma0 = 0.1
+    sigma1 = 0.5
+    #
+    c2 = lambdam * (ffc - ff0) - lambdac * (ffm - ff0)
+    return if c2 >= 0
+        #
+        # Sanity check for negative curvature
+        #
+        lambdap = sigma0 * lambdac
+    else
+        #
+        # It's a convex parabola, so use calculus!
+        #
+        c1 = lambdac * lambdac * (ffm - ff0) - lambdam * lambdam * (ffc - ff0)
+        lambdap = -c1 * 0.5 / c2
+        #
+        lambdaup = sigma1 * lambdac
+        lambdadown = sigma0 * lambdac
+        lambdap = max(lambdadown, min(lambdaup, lambdap))
+    end
+end
+
 function newton_krylov(F, u₀, M::Int = length(u₀); kwargs...)
     F!(res, u) = (res .= F(u); nothing)
     return newton_krylov!(F!, u₀, M; kwargs...)
@@ -195,6 +276,7 @@ function newton_krylov!(
             # ‖F′(u)d + F(u)‖ <= η * ‖F(u)‖ Inexact Newton termination
             kwargs = (; rtol = η, kwargs...)
         end
+        n_res_prior = n_res
 
         # Solve: Jx = -res
         # res is modifyed by J, so we create a copy `-res`
@@ -202,15 +284,26 @@ function newton_krylov!(
         solve!(solver, J, -res; kwargs...)
 
         d = solver.x # Newton direction
-        s = 1        # Newton step TODO: LineSearch
+        # Armijo LineSearch
+        α = 1.0e-4
+        λ₀ = λ = 1.0
+        iarm = -1
+        while true
+            # Update u
+            u .+= λ .* d
+
+            # Update residual and norm
+            F!(res, u) # res = F(u)
+            n_res_armijo = n_res
+            n_res = norm(res)
+            if n_res <= (1 - α * λ) * n_res_prior
+                break # Found success
+            end
 
-        # Update u
-        u .+= s .* d
+            λ = update_lambda(iarm, armfix, λ, λ₀, n_res_prior^2, n_res^2, n_res_armijo^2)
 
-        # Update residual and norm
-        F!(res, u) # res = F(u)
-        n_res_prior = n_res
-        n_res = norm(res)
+            iarm += 1
+        end
 
         if isinf(n_res) || isnan(n_res)
             @error "Inner solver blew up" stats