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 + 2*epsilon*T[x, y], x == xmax] -
  NeumannValue[g + epsilon*q*T[x, y], y == ymax]}, {u,
T}, {x, y} ∈ Ω, {n}];
enter code here