Multiple diffusion reactions ave frequently encountered in the modeling of heterogeneous catalytic reactors. Obtaining an accurate estimate of the yield and selectivity in such reactions is crucial for an optimal design of reactors. Due to the inadequacy of analytical techniques in handling nonuniform catalyst shapes and mixed boundary conditions, numerical techniques are often employed to compute these design parameters. Among other numerical techniques, the boundary element method (BEM) is a superior method to solve linear diffusion reaction problems. The integral nature of the BEM formulation allows for boundary-only discretization of the particle, thus reducing the computer execution time and the data preparation effort. A boundary element algorithm is developed to solve a network of linear diffusion reactions in porous catalyst particles in two dimensions. For this purpose a matrix of fundamental solutions is defined and derived. The developed algorithm is applied to complex reaction networks to obtain the yield of intermediates for nonregular catalyst shapes and nonuniform boundary conditions. The method can be used as a design tool to study particle scale modeling in detail and can be incorporated into an overall reactor model.