The same type of question has been answered here but for single PDE.

The main issue I am facing is how to write the weak form for two PDEs. $ $ \epsilon^2\frac{\partial ^2u}{\partial x^2}+\frac{\partial ^2u}{\partial y^2}=1$ $ $ $ \epsilon^2\frac{\partial ^2T}{\partial x^2}+\frac{\partial ^2T}{\partial y^2}=\epsilon u$ $ $ $ \frac{\partial u}{\partial x}|_{x=0}=0,\,\,\frac{\partial u}{\partial y}|_{y=0}=0,\,u=-2n\frac{\partial u}{\partial y}|_{y=1},\,u=-2\epsilon n\frac{\partial u}{\partial x}|_{x=1}$ $ $ $ \frac{\partial T}{\partial x}|_{x=0}=0,\,\,\frac{\partial T}{\partial y}|_{y=0}=0,\,T=-2\epsilon\frac{\partial T}{\partial y}|_{y=1},\,T=-2\epsilon n\frac{\partial T}{\partial x}|_{x=1}$ $

Any suggestion? or if there is any other way to do it instead of FEM?