With A == [Ax, Ay], etc, and a point P == [Px,Py] for certain u,v, the first step is indeed to subtract A from B, C, D and P, effectively placing A at the origin, i.e, Ax == Ay == 0, just like you do.
Then, you obtain P from u and v as
Px == u*(1 - v)*Bx + (1 - u)*v*Cx + u*v*Dx;
Py == u*(1 - v)*By + (1 - u)*v*Cy + u*v*Dy;
in standard degree 1 Bezier patch notation. The last term includes u*v, which complicates solving for u and v; the solution is indeed nonlinear in Px and Py.
You get
syms Bx By Cx Cy Dx Dy Px Py u v
[solu, solv] = solve([u*(1 - v)*Bx + (1 - u)*v*Cx + u*v*Dx - Px == 0,
u*(1 - v)*By + (1 - u)*v*Cy + u*v*Dy - Py == 0], [u, v])
which you can simplify using the 2D cross product, if you want, like
cross2D = @(AA,BB) AA(1)*BB(2) - AA(2)*BB(1);
syms B C D P
B = [Bx,By];
C = [Cx,Cy];
D = [Dx,Dy];
P = [Px,Py];
% and so on
and then e.g.
syms sigma_2
sigma_2 = 2 * cross2D(C,D-B)
Such cross product are abundant in the solution.
You can then use these results to write a general function uv_from_P(A, B, C, D, P)
