If your system of polynomials is given as a composition of several systems, you can exploit this in HomotopyContinuation.jl.
Composition of two systems
If $f$ and $g$ are your systems, and you want to solve $f \circ g$, then you can do this by
solve(f ∘ g)
For instance, if
$$ f = \begin{bmatrix} ab - 2\\ ac- 1\end{bmatrix}, \quad g = \begin{bmatrix}x + y\\ y + 3\\ x + 2\end{bmatrix},$$
then you solve $f\circ g$ by
julia> @var a b c x y
julia> f = System([a * b - 2, a * c- 1])
julia> g = System([x + y, y + 3, x + 2])
julia> solve(f ∘ g)
Result with 2 solutions
==================================
• 2 non-singular solutions (2 real)
• 0 singular solutions (0 real)
• 4 paths tracked
• random seed: 446964
Composition of several systems
You can also iterate the above process:
julia> @var u v
julia> h = System([u^2 - 1, u + v - 2])
julia> solve(f ∘ g ∘ h)
Result with 4 solutions
==================================
• 4 non-singular solutions (2 real)
• 0 singular solutions (0 real)
• 12 paths tracked
• random seed: 386832
Systems with parameters
You can also track a parameter homotopy in a composite system. If
$$ f = \begin{bmatrix} ab - q\\ ac - p\end{bmatrix}, \quad g = \begin{bmatrix}x + y\\ y + 3\\ x + 2\end{bmatrix},$$
and you want to track solutions from $(p,q) = (1, 2)$ to $(p,q) = (2, 3)$, this is how it works:
julia> res = solve(f ∘ g)
julia> @var p q
julia> f2 = System([a * b - q, a * c - p]; parameters = [p, q])
julia> res2 = solve(f2 ∘ g, solutions(res), p₁=[1, 2], p₀=[2, 3])
Result with 2 solutions
==================================
• 2 non-singular solutions (2 real)
• 0 singular solutions (0 real)
• 2 paths tracked
• random seed: 48821
It does not matter at which level the parameters are.