The numerical investigation of airflow and chemical transport characteristics for a general class of buildings involves identifying values for model parameters, such as effective leakage areas and temperatures, for which a fair amount of uncertainty exists. A Monte Carlo simulation, with parameter values drawn from likely distributions using Latin Hypercube sampling, helps to account for these uncertainties by generating a corresponding distribution of simulated results. However, conducting large numbers of model runs can challenge a simulation program, not only by increasing the need for fast algorithms, but also by proposing specific combinations of parameter values that may define difficult numerical problems. The paper describes several numerical approaches to improving the speed and reliability of the COMIS multizone airflow simulation program. Selecting a broad class of algorithms based on the mathematical properties of the airflow systems (symmetry and positive-definiteness), it evaluates new solution methods for possible inclusion in the COMIS code. In addition, it discusses further changes that will likely appear in future releases of the program.