This study presents a Laplace‐transform‐based boundary element method to model the groundwater flow in a heterogeneous confined finite aquifer with arbitrarily shaped boundaries. The boundary condition can be Dirichlet, Neumann or Robin‐type. The derived solution is analytical since it is obtained through the Green's function method within the domain. However, the numerical approximation is required on the boundaries, which essentially renders it a semi‐analytical solution. The proposed method can provide a general framework to derive solutions for zoned heterogeneous confined aquifers with arbitrarily shaped boundary. The requirement of the boundary element method presented here is that the Green function must exist for a specific PDE equation. In this study, the linear equations for the two‐zone and three‐zone confined aquifers with arbitrarily shaped boundary is established in Laplace space, and the solution can be obtained by using any linear solver. Stehfest inversion algorithm can be used to transform it back into time domain to obtain the transient solution. The presented solution is validated in the two‐zone cases by reducing the arbitrarily shaped boundaries to circular ones and comparing it with the solution in Lin et al. (2016, https://doi.org/10.1016/j.jhydrol.2016.07.028). The effect of boundary shape and well location on dimensionless drawdown in two‐zone aquifers is investigated. Finally the drawdown distribution in three‐zone aquifers with arbitrarily shaped boundary for constant‐rate tests (CRT) and flow rate distribution for constant‐head tests (CHT) are analyzed.