## 1 Introduction  
In this paper, we develop a mathematical model of the filtration flow of a liquid in a piecewise homogeneous porous medium and apply it to the simulation of lymph flow in a lymph node. The classical model of filtration flow relies on [PERSON]'s law, which is valid in a linear porous medium. The numerical methods for solving problems of natural convection in porous media go back to the work of [PERSON] et al. [1], where the finite difference method (FDM) was applied. [PERSON] and [PERSON] [2, 3] used the finite element method (FEM), and [PERSON] and [PERSON] [4] applied the finite volume method (FVM). The finite element and finite volume methods have been widely developed for solving various problems, including those in complex inhomogeneous regions (see, for example, [5, 6, 7, 8]). The traditional application of the filtration flows models is to the problem of fluid flow in different soils. The typical examples are the water intake problem and the oil flow to oil wells. On the other hand, modeling the fluid convection in living organisms and, in particular, the modeling of the lymph flow has its own specifics. There are just a few works that have considered the lymph flow in the lymph node by means of direct numerical simulations. The fluid flow description by [PERSON]'s law in inhomogeneous regions (hydraulic conductivity is not constant) is due to [PERSON] et al. [9]. In [PERSON] et al. [10], a filtration domain was piecewise homogeneous, and the filtration flow was viscous, described by the Darcy-Brinkman law. The absorption of lymph into the blood in both papers was characterized by [PERSON]'s law, which relates the divergence of the velocity field to the fluid pressure. In all these models, the FEM approach was used. If the domain consists of several homogeneous subdomains, then the boundary integral equations method is preferable, and the main goal of this work is to apply this method to the lymph filtration problem. In fact, this method is effective if an integral representation for a solution over domain and subdomains boundaries exists. If the latter is true, then the approach has a number of advantages over the finite difference and finite element methods applied directly in the spatial domain. First, the dimension of the problem being solved is actually reduced. Two-dimensional integral equations (two-dimensional grids) are used instead of three-dimensional grids. Moreover, the differential equations outside the boundary surfaces and conservation laws for the numerical solutions are fulfilled. Finally, the integral representations of unknown functions allow us to calculate derivatives and various functionals without a significant loss of accuracy compared to the accuracy of numerical approximation of the functions themselves. We can single out a number of papers in which the boundary element method was applied to the modeling of flows in porous media governed by Darcy's law or the Darcy-Brinkman law [11]. The solution of a viscous fluid filtration problem governed by the Darcy-Brinkman law using the boundary elements method (BEM) was constructed in [12; 13]. However, these papers only considered two-dimensional problems in a homogeneous domain with the boundary condition set only on the outer boundary. The method of fundamental solutions (MFS) was presented in [14; 15; 16] for two-dimensional problems. This method was similar to the boundary element method, except that the solution was a superposition of partial solutions, the functions of source points located outside the flow area. That allowed avoiding singularity in the boundary condition. In this case, the method of boundary integral equations was applied to flows in one homogeneous region. The application of the boundary integral equations method to the filtration problems in piecewise homogeneous domains was developed in the works of [PERSON] et al., where two-dimensional flows were studied. A systematic presentation of these results can be found in [17].