Question: caculation of an elliptic PDEs system

Hi everyone

I wonder if the following elliptic PDEs system can be caculated by Maple 2018.

The PDEs system is

where

Maple sketched solutions of the above system, but unfortunately main portions (sine portion) of all the solutions seemed to overlap together.

Is it possible to output all the solutions filling the target region (or at least not overlapping) by revising the code?

I am not the sure if the overlap is caused by the boundary conditions(BC) or caused by the elliptic PDEs system??

Particularly, I do not understand how (D[2](x))(xi, 0) = 0 and (D[2](y))(xi, 0) = 0  work in the BCs. Is there any paper or textbook explaining this setup??

Thanks in advance.

The code is as follows (or the code can be downloaded from attatched mw QUESTION_20210423_001.mw):

restart;
with(plots);


alpha04 := (diff(x(xi, eta), eta))^2+(diff(y(xi, eta), eta))^2;

beta04 := (diff(x(xi, eta), xi))*(diff(x(xi, eta), eta))+(diff(y(xi, eta), xi))*(diff(y(xi, eta), eta));

gamma04 := (diff(x(xi, eta), xi))^2+(diff(y(xi, eta), xi))^2;

PDE04 := {alpha04*(diff(x(xi, eta), xi, xi))-2*beta04*(diff(x(xi, eta), xi, eta))+gamma04*(diff(x(xi, eta), eta, eta)) = 0, alpha04*(diff(y(xi, eta), xi, xi))-2*beta04*(diff(y(xi, eta), xi, eta))+gamma04*(diff(y(xi, eta), eta, eta)) = 0};


IBC04 := {x(0, eta) = eta, x(1, eta) = eta, x(xi, 0) = sin((2.0*Pi)*xi), y(0, eta) = eta, y(1, eta) = eta, y(xi, 0) = 100*xi, (D[2](x))(xi, 0) = 0, (D[2](y))(xi, 0) = 0};

pds04 := pdsolve(PDE04, IBC04, numeric, spacestep = 1/50);


xCurve04 := plots[display]([seq(pds04:-plot(x, eta = i), i = 0 .. 10)]);

yCurve04 := plots[display]([seq(pds04:-plot(y, eta = (1/10)*i), i = 0 .. 10)]);