Solve \(\bigtriangledown ^2 T -20 T =-20\) on the unit square \(0\leq x \leq 1,0 \leq y \leq 1\) subject to insulated boundary conditions at \(x=0\) and \(y=0\) and \(T=0\) at \(x=1,y=1\).
In Mathematica 10.0, one can use Finite elements directly. Here is the code to solve the above.
Clear[u,x,y]; r=NDSolveValue[{D[u[x,y],{x,2}]+D[u[x,y],{y,2}]-20*u[x,y]== -200+NeumannValue[0,x==0] +NeumannValue[0,y==0], DirichletCondition[u[x,y]==0,x==1], DirichletCondition[u[x,y]==0,y==1]}, u,{x,0,1},{y,0,1}, Method->{"FiniteElement", "MeshOptions"->{"BoundaryMeshGenerator"->"Continuation"}}];
Plotting the solution
Plot3D[r[x,y],{x,0,1},{y,0,1},Mesh->All]
Now we make a contour plot
ContourPlot[r[x,y],{x,0,1},{y,0,1}]