solve.py

from sympy import Symbol, init_printing
from sympy.solvers import solve
from IPython.display import display

init_printing() 

x1 = Symbol('x1')
x2 = Symbol('x2')
y  = Symbol('y')

p12 = Symbol('p12')
p1y = Symbol('p1y')
p2y = Symbol('p2y')

a = Symbol('a')
b = Symbol('b')
c = Symbol('c')
d = Symbol('d')
e = Symbol('e')
f = Symbol('f')

z0 = Symbol('z0')

sols = solve([
    z0 * a + (1 - z0) * b - x1,
    z0 * c + (1 - z0) * d - x2,
    z0 * e + (1 - z0) * f - y,

    z0 * a * c + (1 - z0) * b * d - p12,
    z0 * a * e + (1 - z0) * b * f - p1y,
    z0 * c * e + (1 - z0) * d * f - p2y
], a, b, c, d, e, f, dict=True)

def compute_sols(x1_, x2_, y_, p12_, p1y_, p2y_, z0_):
    def evaluate(s):
        return s.evalf(subs={x1: x1_, x2: x2_, y: y_, p12: p12_, p1y: p1y_, p2y: p2y_, z0: z0_})
    
    return [dict([(k, evaluate(s)) for k, s in sols[i].items()]) for i in range(2)]

print(compute_sols(0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5))