Solution (source code)

= Solution

Label the sides by $L,R,B,T$, and parametrize the two vertical sides by $s=y\in[-1,1]$, the horizontal sides by $s=x\in[-1,1]$. Their outward <normal derivatives> are $q_L=-u_x(-1,s)$, $q_R=u_x(1,s)$, $q_B=-u_y(s,-1)$ and $q_T=u_y(s,1)$. Let $f_j$ denote the prescribed side values. Assume compatible, sufficiently regular boundary traces; the square's corners have zero arclength measure and do not require separate normal values.

For <sine collocation of square modified Helmholtz global relations>, use <Legendre polynomials> for the known <Dirichlet boundary condition> and a <Fourier sine series> for the unknown <normal derivative>:
$$
f_j(s)\simeq f_j^{(M)}(s)=\sum_{\ell=0}^M d_{j\ell}P_\ell(s),\qquad q_j(s)\simeq q_j^{(N)}(s)=\sum_{m=1}^N c_{jm}\phi_m(s),\qquad \phi_m(s)=\sin[p_m(s+1)],\quad p_m=\frac{m\pi}{2}.
$$
The known coefficients are $d_{j\ell}=(2\ell+1)\int_{-1}^1f_jP_\ell\,ds/2$. The sine functions have $\int_{-1}^1\phi_m\phi_n\,ds=\delta_{mn}$ and form a complete <basis> in $L^2(-1,1)$. Choosing them for the <normal derivative> does not impose zero flux at a corner: the expansion is an $L^2$ representation, and endpoint values are not determined by it. If preserving corner values of the approximated Dirichlet trace is necessary, subtract its endpoint-interpolating line before <polynomial> approximation, and add that line back. All subsequent known-data integrals can alternatively be evaluated with the exact $f_j$.

Define the <entire function> basis transforms
$$
\mathsf P_\ell(\tau)=\int_{-1}^1e^{\tau s}P_\ell(s)\,ds,\qquad \mathsf S_m(\tau)=\int_{-1}^1e^{\tau s}\phi_m(s)\,ds=\frac{p_m[e^{-\tau}-(-1)^me^\tau]}{\tau^2+p_m^2}.
$$
At $\tau=\pm ip_m$ the quotient is evaluated by its removable limit or by the defining integral. For example, $\mathsf P_0=2\sinh\tau/\tau$, and <polynomial> expansion of $P_\ell$ expresses every $\mathsf P_\ell$ in derivatives of $\mathsf P_0$. Put
$$
F_j^{(M)}(\tau)=\sum_{\ell=0}^Md_{j\ell}\mathsf P_\ell(\tau),\qquad Q_j^{(N)}(\tau)=\sum_{m=1}^Nc_{jm}\mathsf S_m(\tau).
$$
For compactness write $A=A(\lambda)$ and $B=B(\lambda)$. Substituting these expansions into the first <global relation for a linear boundary value problem> gives
$$
\boxed{e^{-A}[Q_L^{(N)}(B)+AF_L^{(M)}(B)]+e^A[Q_R^{(N)}(B)-AF_R^{(M)}(B)]+e^{-B}[Q_B^{(N)}(A)+BF_B^{(M)}(A)]+e^B[Q_T^{(N)}(A)-BF_T^{(M)}(A)]\simeq0.}
$$
The signs are fixed by the outward <normal vectors>, rather than by a choice of traversal direction. The second approximate <global relation for a linear boundary value problem> replaces $A$ by $-A$, leaving $B$ unchanged:
$$
\boxed{e^A[Q_L^{(N)}(B)-AF_L^{(M)}(B)]+e^{-A}[Q_R^{(N)}(B)+AF_R^{(M)}(B)]+e^{-B}[Q_B^{(N)}(-A)+BF_B^{(M)}(-A)]+e^B[Q_T^{(N)}(-A)-BF_T^{(M)}(-A)]\simeq0.}
$$
Here $\simeq0$ records omitted boundary-expansion tails. In a finite <spectral method> the selected equations are set equal to zero to solve for the $4N$ unknown real <coefficients> $c_{jm}$. For real data the two families obey the same <complex conjugation> relation as their exact counterparts.