Subsidence due to underground mining may cause damage to surface structures and farmland and change surface and subsurface drainage patterns. During the past few decades, a number of subsidence prediction techniques have been proposed. Currently, the influence function method is the most widely used method for prediction of longwall coal mining subsidence. The influence function method may provide accurate results if a suitable function can be established for a specific site or region based on field observations. The method, however, suffers from generality, and cannot be used in virgin areas. Furthermore, there is a need to develop a model which can simultaneously analyze in-mine stability of mine workings including roof-pillar-floor interaction and associated surface and subsurface subsidence. Since the rock mass surrounding coal seams usually consists of distinct layers/beds and field observations suggest relative movements along bed interfaces, a homogeneous isotropie continuum model cannot effectively describe ground movements due to coal mining. Therefore, in finite element modeling some investigators reduce the values of the Young's modulus to very small, say one percent of the original laboratory value, to obtain acceptable subsidence results, e.g. Bowders et al [1]. Valuable attempts have also been made to model the sliding between beds in finite element modeling [2], but results have been presented only for two-dimensional problems due to large memory and computation time requirements for three-dimensional problems. A three-dimensional numerical model based on a laminated medium used in conjunction with boundary element technique was developed to simulate the ground movements due to coal mining, and results obtained were reasonable [3]. Recently, this model has been further modified to simultaneously analyze both subsidence and in-mine stability of mine workings including roof-pillar-