diff --git a/solve.jl b/solve.jl
index 3c9e5cd..c72eec1 100644
--- a/solve.jl
+++ b/solve.jl
@@ -20,24 +20,26 @@ function solve(F, (G, roots) = start_system(F), maxsteps=10000)
   #  F=homogenize(F)
   H=homotopy(F,G)
   solutions = []
+  steps = 0
 
   @time Threads.@threads for r in roots
     t = 1.0
     step_size = 0.01
     x0 = r
     m = 0
+    steps = 0
 
-    while t > 0 && maxsteps > 0
+    while t > 0 && steps < maxsteps
       x = en_step(H, x0, t, step_size)
       (m, step_size) = adapt_step(x, x0, step_size, m)
       x0 = x
       t -= step_size
-      maxsteps -= 1
+      steps += 1
     end
     push!(solutions, x0)
   end
 
-  return solutions
+  return (solutions, steps)
 end
 
 # Input polynomial system
@@ -46,9 +48,17 @@ F = [x*y - 1, x^2 + y^2 - 4]
 T = [x*y - 1, x^2 + y^2 - 2]
 C = [x^3 - y + 5x^2 - 10, 2x^2 - y - 10]
 
-sF = filter(u -> imag(u[1]) < 0.1 && imag(u[2]) < 0.1, solve(F))
-sT = filter(u -> imag(u[1]) < 0.1 && imag(u[2]) < 0.1, solve(T))
-sC = filter(u -> imag(u[1]) < 0.1 && imag(u[2]) < 0.1, solve(C))
+(sF, sf) = solve(F)
+(sT, st) = solve(T)
+(sC, sc) = solve(C)
+
+println(sf)
+println(st)
+println(sc)
+
+sF = filter(u -> imag(u[1]) < 0.1 && imag(u[2]) < 0.1, sF)
+sT = filter(u -> imag(u[1]) < 0.1 && imag(u[2]) < 0.1, sT)
+sC = filter(u -> imag(u[1]) < 0.1 && imag(u[2]) < 0.1, sC)
 
 # Plotting the system and the real solutions
 ENV["GKSwstype"]="nul"