The indirect boundary element method (IBEM) is formulated for a two-dimensional Stokes flow with moving boundary in the presence of gravity. For simulation of the non-Newtonian fluid flow an iterative scheme for solving a nonlinear system of algebraic equation is used to calculate nonlinear viscous terms. The selected scheme requires values of pseudo forces in each boundary element and each internal collocation point obtained from boundary conditions (for velocities and tractions) and velocity field found at the previous iteration. Boundary element integrals are calculated analytically with a singularity extraction and internal cells integrals are evaluated by using the standard Gaussian quadratures. Internal cells sources are estimated numerically within a finite difference approach to velocity derivatives calculation. The Poiseuille flow of power law fluid in a channel is calculated for a benchmark of the IBEM algorithm. The accuracy and convergence tests are presented for a wide range of power law index (0.2-1.2). In the presence of a free surface, no analytical solution is available. The problem of channel filling with power law fluid is analyzed. The shape and the movement of the free surface are investigated through a special numerical front tracking algorithm based on mesh refinement for the free surface and for the solution domain. Elements of high performance computing are presented. The obtained results coincide with those of other authors. It confirms that the presented algorithm and corresponding software can be applied to study a number of low Reynolds flows, for example in a mould filling technology.