A procedure for the simulation of extraction columns is developed using a program, Thermath, based on the Mathematica (R) computer algebra (CA) software. The procedure requires the simultaneous solution, using the Newton-Raphson method, of mass and energy balances, phase equilibrium equations, and equations needed to match the number of degrees of freedom at each stage. Starting from the equipment model equations, analytical derivatives with respect to the problem unknowns were obtained using CA, and automatically implemented in FORTRAN using Thermath. Similarly, fugacity coefficients and residual enthalpies were derived from the Peng-Robinson and the Soave-Redlich-Kwong equations of state, using either the classical or the Huron-Vidal mixing rules, and also implemented using Thermath. The code was used to simulate the separation of products from an enzymatic reaction using supercritical carbon dioxide. Extraction columns for the systems carbon dioxide-ethanol-water and carbon dioxide-limonene-linalool were also simulated. The results show trends compatible with those that should be expected from supercritical extraction columns. (c) 2004 Elsevier B.V. All rights reserved.