  We present an extension of an algorithm for the classical scalar $p$-Laplace
Dirichlet problem to the vector-valued $p$-Laplacian with mixed boundary
conditions in order to solve problems occurring in shape optimization using a
$p$-harmonic approach. The main advantage of the proposed method is that no
iteration over the order $p$ is required and thus allow the efficient
computation of solutions for higher orders. We show that the required number of
Newton iterations remains polynomial with respect to the number of grid points
and validate the results by numerical experiments considering the deformation
of shapes. Further, we discuss challenges arising when considering the limit of
these problems from an analytical and numerical perspective, especially with
respect to a change of sign in the source term.
