# Solving a coupled PDEs using FEM with Robin conditions

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?

Here is my try but I am not sure whether its correct or not,

a = 1/2; b = 1;
epsilon = a/b;
Ω = Rectangle[{0, 0}, {xmax = 2 b, ymax = 2 a}];
c = -{{epsilon^2, 0}, {0, 1}};

op1 = Div[c.Grad[u[x, y], {x, y}], {x, y}] - 1;
g = 0; q = 1/(2 n);

op2 = Div[c.Grad[T[x, y], {x, y}], {x, y}] - epsilon*u[x, y];

solFEM = ParametricNDSolveValue[{op1 == -NeumannValue[
g + epsilon*q*u[x, y], x == xmax] - NeumannValue[g + q*u[x, y], y == ymax],
op2 == -NeumannValue[g + epsilon*q*T[x, y], x == xmax] -
NeumannValue[g + 2*epsilon*T[x, y], y == ymax]}, {u, T}, {x, y} ∈ Ω, {n}];

Plot3D[solFEM[#][][x, y] & /@ {0.01, 0.05, 0.1} // Evaluate, {x, 0, xmax}, {y, 0, ymax}]

• what is $\epsilon$ is that some small number? – Nasser Nov 24 '17 at 5:17
• @Nasser You can treat it as a parameter. $\epsilon \in (0,1)$. – zhk Nov 24 '17 at 5:23
• It may be useful to add the equations and what else you got so far in Mathematica input form. – user21 Nov 24 '17 at 5:57
• @user21 I have included my try and not sure about it. – zhk Nov 25 '17 at 4:22