I think part of the problem is the K value (usually defined as the ratio of the vapor composition and the liquid composition) is also dependent on bubble and dew point values, so you end up with a "circular" logic/calculation.
Many would say that you find K values (and bubble/dew points) by plugging values into your "black box" (aka "simulator program"), "magic" happens inside of the black box, and out come K values and bubble/dew points. I assume that part of your question is trying to understand the magic that happens inside of these black boxes.
Learning VLE calculations like this can occupy a significant part of a college level course, so I don't think one can expect someone on an internet forum to completely recreate such a course of study. I would strongly suggest that you find a copy of a thermodynamics of phase equilibrium text and study it to understand the details and nuances of performing these calculations.
For introductory purposes, I find that the case using Raoult's law can be very instructive -- at least for introducing the computation steps.
Pi=xi*P0i where Pi is the partial pressure of i, xi is the liquid mole fraction of i, and P0i is the vapor pressure of pure i at T.
P=sum(Pi) -- P is total pressure
yi=Pi/P -- yi is vapor mole fraction of i
Ki=yi/xi
aij=Ki/Kj -- aij is relative volatility of i over j
Since we cannot expect a system like yours to obey Raoult's law, this is only suggested as an introductory computation. From there, you can study the effect of using fugacity and activity coefficients to quantify the non-idealities in the system (for a gamma-phi approach to VLE) or study how to use fugacity coefficients in each phase (phi-phi algorithm) to determine the equilibrium concentrations.