###Problem 13.33 Stephen Lincoln, Amanda Card, John Dailey, Amanda Masotta###
import math
from scipy import optimize


def func(x):
    total = (1 - 2 * x[0])
    yh2 = (.75 - 2 * x[0] - x[1]) / total
    yco = (.15 - x[0] + x[1]) / total
    yco2 = (.05 - x[1]) / total
    ych3oh = (x[0]) / total
    yh2o = x[1] / total

    k1 = 6.749 * 10 ** -4
    k2 = 0.01726

    out = [(ych3oh / ((yh2) ** 2 * yco)) - (100 / 1) ** 2 * k1]   ###Equation 1, set eq to zero
    out.append((yco * yh2o / (yh2 * yco2) - k2))                  ###Equation 2, set eq to zero (P/Po)^0 = 1
    return out

extents = optimize.fsolve(func, [.1, .1])  ###initial guess .1 for E1, E2
print("Extents: ", extents)