Can you try to run this code?
using SymPy
@syms x y z A B C D E F a b c d e f g
eq1 = x * (x - 1) * (D + E * y * z - E * y - E * z - A - B * y * z + B * y +
B * z - f * y - g * z - c * y)
eq2 = y * (y - 1) * (F + f * x - f + c * x)
eq3 = -z * (z - 1) * (-b + d + e)
solve([eq1, eq2, eq3], [x, y, z], domain=SymPy.sympy.S.Reals)
I try but I don’t know if the battery lasts long enough 
Why does it takes so long time? I have run for one day, but it still not finishes.
I have no idea.
This is my first time using these packages.
I believe that the use of all these parameters commits (unnecessarily) the algorithm that “rearranges” the equations, for a long time.
If you do it by hand (even the ncomplete system can be solced by hand) you save most of the effort.
I wonder if with some suitable suggestions (via some kwarg) sympy is not able to do by itself the simplifications that I did.
Can you help me correct the codes? The solve can not finish in my computer, I can not get the results. 
you could define A, B, C etc as a function of the original parameters.
Solve the systemq of equations and then do the inverse substitutions.
Since eq1, eq2 are somewhat decoupled from eq3:
julia> solve([eq1,eq2], [x,y], domain=SymPy.sympy.S.Reals)
5-element Vector{Tuple{Sym, Sym}}:
(0, 0)
(0, 1)
(1, 0)
(1, 1)
((-F + f)/(c + f), (-A + B*z + D - E*z - g*z)/(B*z - B - E*z + E + c + f))
terminates within reasonable time. Adding the eq3 consequences isn’t hard.
How to do the inverse substitutions?
using SymPy
@syms x y z e₁ e₂ c₁ c₂ c₃ c₅ s₁ s₂ s₃ g₁ g₂ A B C D E F G
eq1 = x*(x - 1)*(D - y*( A + z * B) - z * C)
eq2 = y*(y - 1)*( E + F * x)
eq3 = -z*(z - 1)*G
solve([eq1, eq2, eq3], [x, y, z])
10-element Vector{Tuple{Sym, Sym, Sym}}:
(0, 0, 0)
(0, 0, 1)
(0, 1, 0)
(0, 1, 1)
(1, 0, 0)
(1, 0, 1)
(1, 1, 0)
(1, 1, 1)
(-E/F, D/A, 0)
(-E/F, (-C + D)/(A + B), 1)
A => c₂ + g₁ + s₁ - e₂
B => e₂ - c₂
C => c₂ + g₂ - e₂
D => c₁ - e₁
E => c₃ - g₁
F => g₁ + s₁
G => -c₅ + s₂ + s₃
I try to use subs, but it gets errors.
using SymPy
@syms x y z e₁ e₂ c₁ c₂ c₃ c₅ s₁ s₂ s₃ g₁ g₂ A B C D E F G
eq1 = x * (x - 1) * (D - y * (A + z * B) - z * C)
eq2 = y * (y - 1) * (E + F * x)
eq3 = -z * (z - 1) * G
answ = solve([eq1, eq2, eq3], [x, y, z])
subs(answ, A => c₂ + g₁ + s₁ - e₂, B => e₂ - c₂, C => c₂ + g₂ - e₂, D => c₁ - e₁, E => c₃ - g₁, F => g₁ + s₁)
ERROR: MethodError: no method matching subs(::Vector{Tuple{Sym, Sym, Sym}}, ::Pair{Sym, Sym}, ::Pair{Sym, Sym}, ::Pair{Sym, Sym}, ::Pair{Sym, Sym}, ::Pair{Sym, Sym}, ::Pair{Sym, Sym})
Closest candidates are:
subs(::T, ::Pair...; kwargs...) where T<:SymPy.SymbolicObject at ~/.julia/packages/SymPy/7cg6m/src/utils.jl:126
subs(::Tuple{T, N}, ::Any...; kwargs...) where {T<:SymPy.SymbolicObject, N} at ~/.julia/packages/SymPy/7cg6m/src/utils.jl:127
subs(::Number, ::Any...; kwargs...) at ~/.julia/packages/SymPy/7cg6m/src/utils.jl:128
...
Stacktrace:
[1] top-level scope
@ ~/Desktop/sfby.jl:44
after some trial and error, this seems to work.
But having absolutely no experience with these packages, I believe there are “easier” solutions
julia> sol=solve([eq1, eq2, eq3], [x, y, z])
10-element Vector{Tuple{Sym, Sym, Sym}}:
(0, 0, 0)
(0, 0, 1)
(0, 1, 0)
(0, 1, 1)
(1, 0, 0)
(1, 0, 1)
(1, 1, 0)
(1, 1, 1)
(-E/F, D/A, 0)
(-E/F, (-C + D)/(A + B), 1)
julia> repl=(A => c₂ + g₁ + s₁ - e₂,
B => e₂ - c₂,
C => c₂ + g₂ - e₂,
D => c₁ - e₁,
E => c₃ - g₁,
F => g₁ + s₁,
G => -c₅ + s₂ + s₃)
(A => c₂ - e₂ + g₁ + s₁, B => -c₂ + e₂, C => c₂ - e₂ + g₂, D => c₁ - e₁, E => c₃ - g₁, F => g₁ + s₁, G => -c₅ + s₂ + s₃)
julia> map(s->foldl((sc,c)->subs.(sc,c),repl, init=s), sol)
10-element Vector{Tuple{Sym, Sym, Sym}}:
(0, 0, 0)
(0, 0, 1)
(0, 1, 0)
(0, 1, 1)
(1, 0, 0)
(1, 0, 1)
(1, 1, 0)
(1, 1, 1)
(-(c₃ - g₁)/(g₁ + s₁), (c₁ - e₁)/(c₂ - e₂ + g₁ + s₁), 0)
(-(c₃ - g₁)/(g₁ + s₁), (c₁ - c₂ - e₁ + e₂ - g₂)/(g₁ + s₁), 1)
That works. The issue is the solution is a vector of tuples. Substitution must be applied to element of each tuple, as it is currently defined. Maybe this pattern seems more direct:
[subs.(s, repl...) for s ∈ sol]
As for speed, I found that by making the variables real (@syms x::real ... or x, y,... = symbols("x, y...", real=true)) the solving went fairly quickly.
Although I can do some substitutions and then finish the solve, but I think this is not a good method. I hope it can directly solve the problem.
You ask in the wrong place. SymPy is not a Julia library, but a Python library. You got very good answers how to work around the issues with it, but if you suggest to improve it please create a ticket at Issues · sympy/sympy · GitHub .