Solution to combined confined and unconfined flows in multiple aquifer systems has been a challenging issue for groundwater flow modelers. In this study, steady-state solutions to four different problems are presented using analytical and numerical methods. For analytical solution, the concept of discharge potential is used. In this method, the governing equation becomes the same for both confined and unconfined flows. Moreover, discharge potential satisfies Laplace's Equation and hence the potential theory becomes applicable which makes easier to reach solutions to complex problems. The problem consisting of an aquifer with a discontinuous aquiclude was solved by using comprehensive potential which gives the solution in an entirely unconfined aquifer with equivalent boundary conditions. This facilitates the solutions to confined and unconfined regions above and below the aquiclude. For the numerical solutions, two different versions of MODFLOW, namely, MODFLOW-2005 and MODFLOW-NWT were used. Initially, two different resolution schemes were tested. One of the schemes yielded good predictions of heads for different types of input in which the flow direction and aquifer type varied.