How do I estimate the number of theoretical stages using the McCabe-Thiele method for a non-ideal system?

The classic McCabe-Thiele method (operating lines plus equilibrium curve, stepping off stages between them) technically assumes constant molar overflow, which strictly only holds for near-ideal systems with similar molar heats of vaporization. For a genuinely non-ideal system, you have a couple of options depending on how far off ideal you are.

If the non-ideality is moderate, you can often still use McCabe-Thiele but generate the equilibrium curve (x-y diagram) from actual VLE data or an activity coefficient model (NRTL, UNIQUAC, Wilson) rather than assuming Raoult's law that alone captures a lot of the real behavior, including azeotropes, which will show up clearly as the equilibrium curve crossing the 45° line. The graphical stepping-off procedure itself still works fine as long as your equilibrium curve is accurate.

If the system also has significantly different heats of vaporization between components (so constant molar overflow genuinely fails), you need to move to an energy-balance-based method - either the Ponchon-Savarit method (which uses an enthalpy-composition diagram and accounts for varying liquid and vapor flows tray to tray), or more practically today, just build it in Aspen Plus or HYSYS with a rigorous tray-by-tray simulation (RadFrac or the equivalent) using a proper activity coefficient or equation-of-state model. Honestly, for anything beyond a classroom problem or first-pass estimate, rigorous simulation has mostly replaced hand methods for non-ideal systems - use McCabe-Thiele for intuition and quick checks, but verify with simulation before you commit to a design.

You can learn the Mccabe-Thiele diagram here - https://youtu.be/gw2rtssRyvo?si=1AZ5S2MWQm0f5jIY