This question is in some sense a continuation of this question, though more focused in its scope.
Let $(M,g)$ be an $n$-dimensional, connected, compact Riemannian manifold with boundary. Assume we are given an immersion $f:M \to \mathbb{R}^n$. (i.e $df$ is invertible at every point $p \in M$, note that I assume $n$ is the dimension of $M$).
Let $\omega:(M,g) \to (\mathbb{R}^n,e)$ be the harmonic function corresponding to the Dirichlet problem imposed by $f$, i.e $\omega|_{\partial M}=f|_{\partial M}$
Is it true that $\omega$ must be an immersion? Does anything change if we assume $f$ is (globally) injective?
Note that since $f$ is an immersion, it is locally injective, hence $f|_{\partial M}$ is locally injective. (In particular naive "counter-examples" like taking $\,f|_{\partial M}$ to be constant so $\omega$ is constant do not work). In fact we have $\text{rank}(d\omega_p)\ge \text{rank}\big(d(\omega|_{\partial M})_p\big)= \text{rank}\big(d(f|_{\partial M})_p\big)=n-1$ for every $p \in \partial M$
Here's a counterexample on the unit disc. Let $f : \mathbb D^2 \to \mathbb R^2$ be defined by
$$ f(x,y) = (x-2y^2,y). $$
This has differential $$\left(\begin{matrix}1 & -4y \\ 0 & 1\end{matrix}\right)$$
and thus is an immersion. It's also injective.
The solution to the Dirichlet problem in this case is $$\omega(x,y) = (x^2 - y^2 + x - 1,y);$$ you can check that this is harmonic and agrees with $f$ on the boundary. However, its differential $$\left(\begin{matrix}1+2x & -2y \\ 0 & 1\end{matrix}\right)$$ is degenerate along the line $x = -1/2$, so $\omega$ is not an immersion.
It seems that you at least need to assume that $f(M)$ is convex - in the two-dimensional case this is enough (see https://eudml.org/doc/152349), but I'm not sure whether such a result holds in higher dimensions.