I think the difficulty with the subdifferential approach you describe lies in finding critical points. For instance, with your example, solving the equations you get by setting partials to $0$ yields only $x_1 = | x_2 |$ and $\lambda = -1$. This set of equations is too undetermined to yield a solution. As you point out, if you were given a solution $x_1 = x_2 = 0$, $\lambda = -1$, you could verify that it is a critical point with the subdifferential approach, but that's quite different from deriving $x_1 = x_2 = 0$, $\lambda = -1$ as a critical point.
Remember, though, that the definition of a critical point of a function includes points at which the derivative fails to exist. So when you're constructing your set of critical points using Lagrange multipliers, look not only for points at which the partials are zero but also for points at which a partial does not exist. For instance, it's clear in your example that the only point at which $\frac{\partial L}{\partial x_2}$ fails to exist is $x_2 = 0$. (You also get $\lambda = 0$ from $\frac{\partial L}{\partial x_2}= 0$, but this is inconsistent with the $\lambda = -1$ from $\frac{\partial L}{\partial x_1}= 0$.) Including that with your equations from $\frac{\partial L}{\partial x_1} = 0$ and $\frac{\partial L}{\partial \lambda} = 0$ quickly gets you the solution $x_1 = x_2 = 0$, $\lambda = -1$.
In fact, finding points where a partial does not exist is generally not any harder (and is sometimes easier) than finding points where a partial is zero. There aren't many ways for a point in the domain of a continuous function to have a partial that fails to exist. (Besides the absolute value situation, you could have division by zero with the partial from an expression like $x^p$, with $0 < p < 1$, in the function. Perhaps there are some others I'm forgetting - maybe others reading this can supply them. There are also some pathological cases, like the continuous but nowhere differentiable Weierstrass function, but those don't generally show up in practice.)
And, of course, if you have more than just a few variables, you don't want to use Lagrange multipliers anyway: Solving a large number of nonlinear equations is just too difficult. In that case, you will probably want an iterative nonlinear optimization method.
Note that
$$T(x,y)=11-(x-1)^2-4y^2\ .$$
This shows that $T(x,y)\leq T(1,0)=11$ for all $(x,y)\in D$, and the maximum value $11$ is taken only at the point $(1,0)$.
It is obvious that $(1,0)$ is the only stationary point of $T$ in the interior of $D$. Therefore the minimum of $T$ has to occur on the boundary $\partial D$. In order to find it we parametrize $\partial D$ by means of
$$\partial D:\quad\phi\mapsto(2\cos\phi,2\sin\phi)\qquad(\phi\in{\mathbb R})$$
and consider the pullback
$$f(\phi):=T(2\cos\phi,2\sin\phi)=\ldots=12\cos^2\phi+4\cos\phi-6\ .$$
The derivative
$$f'(\phi)=-(24\cos\phi+4)\sin\phi$$
vanishes when $\cos\phi=-{1\over6}$ or $\sin\phi=0$. This leads to $4$ conditionally stationary points of $T$ on $\partial D$, namely
$$\left(-{1\over3},{1\over3}\sqrt{35}\right),\quad \left(-{1\over3},-{1\over3}\sqrt{35}\right),\quad (2,0),\quad(-2,0)\ .\tag{1}$$
When $\cos\phi=-{1\over6}$ then $f(\phi)=-{19\over3}$, and when $\sin\phi=0$ then $f(\phi)=6\pm 4>0$. It follows that the minimum of $T$ on $D$ is $-{19\over3}$, and it is taken at the first two points listed in $(1)$.
Best Answer
The extremal temperatures are taken either in the interior $D$ or on the boundary $\partial D$ of the unit disk. Writing $$T(x,y)=2x^2+\left(y-{1\over2}\right)^2-{1\over4}$$ we see that the graph of $T$ is a paraboloid with its apex at $P_1=\bigl(0,{1\over2}\bigr)\in D$, and there are no other stationary points of $T$ in $D$. Take note of $T(P_1)=-{1\over4}$.
For $\partial D$ we use the parametrization $\phi\mapsto(\cos\phi,\sin\phi)$ and then have to consider the pullback $$\hat T(\phi)=2\cos^2\phi+\sin^2\phi-\sin\phi={9\over4}-\left(\sin\phi+{1\over2}\right)^2\ .$$ The function $\hat T$ is maximal (namely $={9\over4}$) at the points where $\sin\phi=-{1\over2}$.
This leads to the candidates $P_2=\bigl(\cos{7\pi\over6},\sin{7\pi\over6}\bigr)$ and $P_3=\bigl(\cos{11\pi\over6},\sin{11\pi\over6}\bigr)$. Take note of $T(P_2)=T(P_3)={9\over4}$.
The function $\hat T$ is minimal at the points where $\sin\phi=1$. This leads to the candidate $P_4=(0,1)$ with $T(P_4)=\hat T\bigl({\pi\over2}\bigr)=0$.
Summing it up we can say that it is coldest (namely $-{1\over4}$) at $P_1\in D$ and hottest (namely ${9\over4}$) at the points $P_2$, $P_3\in\partial D$.