solve_many.jl

using HomotopyContinuation, LinearAlgebra

@polyvar x y z p[1:3]

F = [
    x + 3 + 2y + 2 * y^2 - p[1],
    (x - 2 + 5y) * z + 4 - p[2] * z,
    (x + 2 + 4y) * z + 5 - p[3] * z,
]

# Generate generic parameters by sampling complex numbers from the normal distribution
p₀ = randn(ComplexF64, 3)
# Compute all solutions for F_p₀
S = solutions(solve(subs(F, p => p₀)))
# create some fake paramter values
params = [randn(3) for _ = 1:1000]
# create a `Solver` to reuse for the path tracking
F_solver = solver(F; parameters = p, generic_parameters = p₀)
# solve the system F_p for all parmaters p in params
params_solutions = map(params) do p
    solutions(solve(F_solver, S_p₀; target_parameters = p))
end