How can I make NDSolve aware of vectorial nature of an ODE?

Question

I am trying to solve the equilibrium equation $\text{DIV } \mathbf P(\textbf{u}(x,y)) = \mathbf{f}$, where $\mathbf P$ is the stress tensor defined by

$$ \mathbf P = \mathbf F + (1-\det\mathbf F)^2\quad\text{with}\quad \mathbf F = \text{GRAD } \mathbf u + \mathbf I. $$

I implemented the equations above in Mathematica as follows

ns[x,y] = {nx[x,y], ny[x,y]};
F = Inactive[Grad][ns[x,y], {x,y}] + IdentityMatrix[2];
J = Det[F] //FullSimplify;
μ = 10;
λ = 5;
P = λ F + μ (1 - J)^2
solInt = y * x^2;
equation = Inactive[Div][P, {x, y}] == Inactive[Grad][solInt, {x, y}] ;
boundaryConds1 = {ns[x, y] == {3, 0}};
boundaryConds2 = {ns[x, y] == {0, 0}};
boundaryConds3 = {ns[x, y] == {0, 0}};
Ω = Rectangle[{0, 0}, {1, 1}]
un = NDSolveValue[{equation,
   { 
      DirichletCondition[boundaryConds1, x == 1],
      DirichletCondition[boundaryConds2, x == 0],
      DirichletCondition[boundaryConds3, y == 0]
    }}, ns[x, y], {x, y} ∈ Ω]

This gives me the error:

"There are more dependent variables, {nx[x, y], ny[x, y]}, than equations, so the system is underdetermined."

I suspect that the culprit is the fact that Mathematica does not understand that the equation, once evaluated, will yield a vector with two entries.

Any way to fix the issue?