A numerical procedure for analyzing two-phase flow and heat transfer in porous media was presented. Formulation was based on the two-phase mixture model, originally developed by Wang and Beckermann, in conjunction with a finite element method (FEM). An extended Forchheimer - Darcy law was employed to accommodate the procedure to a wide range of Reynolds number. The conservation equations were discretized by FEM on quadrilateral elements using an equal-order interpolation. To validate the proposed algorithm, boiling in a square, two-dimensional porous enclosure heated from below was simulated. Numerical results agreed well with those by Wang and Beckermann and experimental data by Sondergeld and Turcotte. Then free convection boiling from a horizontal cylinder embedded in porous medium and heated at uniform heat flux was analyzed. Numerical results include the variations in isotherm and iso-liquid saturation, and the flow fields as a function of heat flux and temperature of top surface. The present algorithm is adaptable to a variety of geometric configurations, thermal and hydrodynamic boundary conditions, and a wide range of flow velocity, thus applicable in tackling any types of problems of practical importance.