diff --git a/random_poly.jl b/random_poly.jl
new file mode 100644
index 0000000..f3b7bad
--- /dev/null
+++ b/random_poly.jl
@@ -0,0 +1,26 @@
+module RandomPoly
+  export random_system
+
+  using TypedPolynomials
+  using Random
+  using Distributions
+
+  # Random polynomial of degree n in m variables
+  function random_poly(n, m)
+    x = [TypedPolynomials.Variable{Symbol("x[$i]")}() for i in 1:m]
+
+    monomial_powers=vcat(collect(Iterators.product([0:n for _ in 1:m]...))...)
+
+    return sum(map(i -> rand(Uniform(-10,10)) * prod(x.^i), monomial_powers))
+  end
+
+  # Generate a system of m random polynomials in m variables of degree d_i
+  function random_system(m, max_degree)
+    d = rand(1:max_degree, m)
+    println("generating system")
+    random_polys = [random_poly(d[i], m) for i in 1:m]
+    println("done generating system")
+
+    return random_polys
+  end
+end
diff --git a/report/report.pdf b/report/report.pdf
index 98dfb88..b61f17b 100644
Binary files a/report/report.pdf and b/report/report.pdf differ
diff --git a/solve.jl b/solve.jl
index f11c4be..c84cee8 100644
--- a/solve.jl
+++ b/solve.jl
@@ -2,27 +2,28 @@
 using TypedPolynomials
 
 # Local dependencies
+include("random_poly.jl")
 include("start-system.jl")
 include("homotopy.jl")
-#  include("homogenize.jl")
 include("euler-newton.jl")
 include("adapt-step.jl")
 include("plot.jl")
+using .RandomPoly
 using .StartSystem
 using .Homotopy
-#  using .Homogenize
 using .EulerNewton
 using .AdaptStep
 using .Plot
 
 # Main homotopy continuation loop
 function solve(F, (G, roots) = start_system(F), maxsteps = 1000)
-  #  F=homogenize(F)
   H=homotopy(F,G)
   solutions = []
   step_array = []
 
   Threads.@threads for r in roots
+  #  for r in roots
+    println("New root")
     t = 1.0
     step_size = 0.01
     x0 = r
@@ -43,30 +44,37 @@ function solve(F, (G, roots) = start_system(F), maxsteps = 1000)
 end
 
 # Input polynomial system
-@polyvar x y
-C = [x^3 - y + 5x^2 - 10, 2x^2 - y - 10]
-Q = [x^2 + 2y, y - 3x^3]
-F = [x*y - 1, x^2 + y^2 - 4]
-T = [x*y - 1, x^2 + y^2 - 2]
+dimension = 4
+max_degree = 3
+R = random_system(dimension, max_degree)
+#  @polyvar x y
+#  C = [x^3 - y + 5x^2 - 10, 2x^2 - y - 10]
+#  Q = [x^2 + 2y, y - 3x^3]
+#  F = [x*y - 1, x^2 + y^2 - 4]
+#  T = [x*y - 1, x^2 + y^2 - 2]
 
-(sC, stepsC) = solve(C)
-(sQ, stepsQ) = solve(Q)
-(sF, stepsF) = solve(F)
-(sT, stepsT) = solve(T)
+(sR, stepsR) = solve(R)
+#  (sC, stepsC) = solve(C)
+#  (sQ, stepsQ) = solve(Q)
+#  (sF, stepsF) = solve(F)
+#  (sT, stepsT) = solve(T)
 
-println("C: ", stepsC)
-println("Q: ", stepsQ)
-println("F: ", stepsF)
-println("T: ", stepsT)
+println("R: ", stepsR)
+println("solutions:" sR)
 
-sC = filter(u -> imag(u[1]) < 0.1 && imag(u[2]) < 0.1, sC)
-sQ = filter(u -> imag(u[1]) < 0.1 && imag(u[2]) < 0.1, sQ)
-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)
+#  println("C: ", stepsC)
+#  println("Q: ", stepsQ)
+#  println("F: ", stepsF)
+#  println("T: ", stepsT)
+
+#  sC = filter(u -> imag(u[1]) < 0.1 && imag(u[2]) < 0.1, sC)
+#  sQ = filter(u -> imag(u[1]) < 0.1 && imag(u[2]) < 0.1, sQ)
+#  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)
 
 # Plotting the system and the real solutions
-ENV["GKSwstype"]="nul"
-plot_real(sC, C, 6, 12, "1")
-plot_real(sQ, Q, 2, 2, "2")
-plot_real(sF, F, 4, 4, "3")
-plot_real(sT, T, 4, 4, "4")
+#  ENV["GKSwstype"]="nul"
+#  plot_real(sC, C, 6, 12, "1")
+#  plot_real(sQ, Q, 2, 2, "2")
+#  plot_real(sF, F, 4, 4, "3")
+#  plot_real(sT, T, 4, 4, "4")