I doubt that the original poster knows anything about measure theory, much less Hausdorff measure.

Yes, Crank-Nicholson is unconditionally stable. The naive FTCS (forward differences in time, centered differences in space) scheme isn't stable unless you're careful about picking $\Delta x$ and $\Delta t$. If you just want to solve an equation, then by all means, use Crank-Nicholson, but if you want to understand what's going on, then the FTCS scheme is a good place to start.

You're right of course that from a physical point of view this isn't too surprising, since heat diffusion is about random vibration at the molecular scale. However, if you start with the mathematical problem $u_{t}=u_{xx}$, it's quite surprising that the solution should have any connection to probability theory. Most of my students (mathematics majors in a senior modeling course) have limited physics backgrounds and are surprised by this.

My answer may not have been clear about one important point- I'm talking about the time dependent heat equation (which is a parabolic BVP.) The elliptic version of the problem in which you find the steady state solution is easier but not as interesting. Also, in my experience of teaching this to students at your level, it's been considerably easier to start with the problem in one space dimension and time rather than starting in two dimensions.

I'd argue that this wasn't so much a matter of intuition as knowing some tricks that are frequently useful in convex optimization. I was very familiar with the idea of using $\mbox{vec}()$ to convert the Frobenius norm of a matrix into the 2-norm of a vector and with the idea of writing the the matrix triple product with a diagonal matrix in the middle as a sum of outer products. If you'd like to learn more of this, I'd strongly encourage you to read the textbook "Convex Optimization" by Vandenberghe and Boyd.

Yep, you've got to give up on $C\geq 0$.

The motivation is to turn this problem into a semidefinite programming problem. This is a standard construction in SDP.

Just to be clear, you want to find a matrix $M$ (linearly dependent on $B$), so that $M \ge 0$ iff $A-BC^{-1}B^{T} \leq 0$, right?

No, I didn't say that. You've made a sign change error in what you wrote: $A-BC^{-1}B^{T} \leq 0$ iff $-A+BC^{-1}B^{T} \geq 0$. The signs on both terms have to change.