This thesis is focused on developing a numerical method for solving a problem where fluid flow governed by Stokes equations is coupled with flow through a porous medium which is governed by Darcy's equations. For closed smooth boundaries, flow in each region is represented as boundary integrals with densities to be determined to satisfy the boundary conditions. With the solution in this form the dimension of the problem reduces but the integrands are singular. The proposed numerical method is based on regularizing the integrands and discretizing the integrals by a quadrature. The numerical error is discussed and some convergence results are shown for Stokes flow in 3D. A way of reducing the error is discussed. Next, the fluid quantities in two dimensions are reformulated as single and double layer potentials and a solution method with higher accuracy is proposed. It is based on regularizing the kernels and subtracting the highest regularization error term. The method is applied to several test cases of Stokes and Darcy flows, and an increase of the convergence rate from first to second is observed. In the coupled case, the two solutions are coupled through appropriate interface conditions, which in this case do not require the velocity to be continuous. An example of a coupled problem is presented A fast summation technique for impulse methods in three dimensions is presented. Impulse methods provide a representation of Euler flows in terms of impulse variables. The fundamental solution of Darcy flow is a potential dipole, thus this is a robust method for computing Darcy velocity for dense grids