This question is a follow-on of this answer of mine. This answer provides a numerical solution ; the initial question was about possible closed form expressions.
Q1 : Is there a reference in the existing literature about these orthogonal curves ?
Q2 : These orthogonal curves look symmetrical with respect to second diagonal line with equation $x+y=1$. But is it the case ?
Q3 : Let us consider one of the orthogonal functions $(x(t),y(t))$. Let $z(t)$ be any or $x(t)$, $y(t)$ or $z(t)=x(t)+i y(t)$. This function verifies (by converting the first order system into a second-order ODE) :
$$z''(t)+\frac{2}{\sin(2t)}z'(t)+z(t)=0 \tag{*}$$
(with appropriate initial conditions) or the equivalent equation written under the Sturm-Liouville form :
$$\frac{d}{dt}(\tan(t) \ z'(t))+\tan(t) \ z(t)=0 \tag{**}$$
but with eigenvalue $\lambda = 0$ (?). Is it possible to take advantage of these forms (*) or (**) in order to have a different approach to this issue ?
Remark : If the initial points $I_k$ (notations of my text) are replaced by their symmetrical points on the "Eastern" side of square $S$
(symmetry with respect to diagonal line with equation $x+y=1$) and if we replace instruction yp=[Y/ta;X*ta] by yp=[-Y/ta;-X*ta], we get (numericaly speaking) the same curves.