save some differentiate calls, fix newton residual check

euler-newton.jl

@@ -11,19 +11,21 @@ module EulerNewton
     vars = variables(H(t))
     # Jacobian of H evaluated at (x,t)
     JH = [jh(vars=>x) for jh in differentiate(H(t), vars)]
-    Δx = JH \ -[gg(vars=>x) for gg in H(1)-H(0)] # ∂H/∂t is the same as γG-F=H(1)-H(0) for our choice of homotopy
+    # ∂H/∂t is the same as γG-F=H(1)-H(0) for our choice of homotopy
-    xp = x .+ Δx * step_size
+    Δx = JH \ -[gg(vars=>x) for gg in H(1)-H(0)]
+    xh = x + Δx * step_size
 
     # Corrector step
+    JHh=differentiate(H(t+step_size), vars)
     for _ in 1:10
-      JH = [jh(vars=>xp) for jh in differentiate(H(t+step_size), vars)]
+      JH = [jh(vars=>xh) for jh in JHh]
-      Δx = JH \ -[h(vars=>xp) for h in H(t+step_size)]
+      Δx = JH \ -[h(vars=>xh) for h in H(t+step_size)]
-      xp = xp .+ Δx
+      xh = xh + Δx
-      if LinearAlgebra.norm(Δx) < 1e-6
+      if LinearAlgebra.norm([h(vars=>xh) for h in H(t+step_size)]) < 1e-8
         break
       end
     end
 
-    return xp
+    return xh
   end
 end