Coupling the adjacent zones for seepage analysis in porous media needs compatibility and equilibrium equations (equality of potential on coinciding nodes and conservation of flowing mass between zones, respectively). When stretched coordinate transformation is applied to the anisotropic zones, the Dirichlet boundary conditions remain unchanged, but the Neumann boundary condition should also be transformed. Similarly in a zoned problem, for the interface between zones, compatibility equations remain unchanged during the transformation while the equilibrium equations should be transformed. In this paper, transformed Neumann boundary conditions and equilibrium equations for the interface of neighbor anisotropic zones for seepage problems have been developed in three dimensions. A computer program for seepage analysis of zoned anisotropic media based on the Boundary Element Method is developed. The code is used to solve several examples with isotropic and anisotropic zones. Some examples are also solved by finite element method for verification. Illustrated results show the ability and accuracy of the mathematical and the numerical model for solving different types of applied three-dimensional seepage problems that arise in engineering practice.