A mathematical model has been developed to predict the steady state performance of a continuous bulk styrene polymerization process with catalytic initiation for solid polystyrene. The polymerization section contains one boiling CSTR followed by multiple linear-flow reactors. The devolatilization section consists of two polymer pre-heaters and two high-solids flashes. The polymer moment equations were solved simultaneously with the reactor modeling equations. The non-linear algebraic equations were solved by a Newton-Raphson iteration technique to give the steady-state styrene monomer weight fraction in a CSTR. The coupled, non-linear ordinary differential equations were numerically integrated using a single-step, 4th-order Runge-Kutta technique, followed by a multi-step Adams-Moulton technique. The resulting computer simulation model is capable of evaluating how the production rate and product quality are affected by feed composition, temperature, initiator type, initiator concentration, and residence time. Several case studies were given for commercially important crystal-clear and impact-resistant resins. A binary initiation system gives a good balance of monomer conversion, polymer molecular weights, and rubber grafting compared to a single initiation system. The styrene dimer/trimer occur in low concentrations but can be substantially reduced with a low temperature initiator. The ideal mean residence time is approximately one minute or less in a shell-and-tube devolatilization pre-heater. Low flash chamber vacuum is more effective than high polystyrene melt temperature to reduce the volatile content of the final product. The water injected to the low volatile melt shows promising improvement in the second-stage polystyrene devolatilization.