We consider shape optimization of Stokes flow in channels where the objective is to design the lateral walls of the channel in such a way that a desired velocity profile is achieved. This amounts to the solution of a PDE constrained optimization problem with the state equation given by the Stokes system and the design variables being the control points of a Bézier curve representation of the lateral walls subject to bilateral constraints. Using a finite element discretization of the problem by Taylor–Hood elements, the shape optimization problem is solved numerically by a path-following primal-dual interior-point method applied to the parameter dependent nonlinear system representing the optimality conditions. The method is an all-at-once approach featuring an adaptive choice of the continuation parameter, inexact Newton solves by means of right-transforming iterations, and a monotonicity test for convergence monitoring. The performance of the adaptive continuation process is illustrated by several numerical examples.