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

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-detmathbf F)^2quadtext{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?