(Subscript[u, t])=k Δu
Radially Symmetric on 
0<r<rMax and 0<θ<2π

We now have three variables θ,r, and t.  As of now we only know how to do two variables.  We have a choice:
r an,

 is a time varying radial (constant wrt angle θ) temperature distribution;
r and θ is an equilibrium (constant wrt time) temperature distribution;
θ and t is a time varying angular (constant wrt radius) temperature distribution.

(Subscript[u, t])=k Δu with u(Subscript[|, ∂ω])=0
0<r<rMax and 0<θ<2π

We have three variables θ,r, and t.  This is our first three variable problem.  The separation trick still works it just gets a little messier. 

As before assume that u=T[t] ϕ[r,θ] and plug in to the PDE 
	({{T'[t]ϕ[r,θ], =, k T[t] Δϕ[r,θ]}})
and separate to get
	({{((T'[t])/(k T[t])), =, ((Δϕ[r,θ])/(ϕ[r,θ]))}})=-λ
The time equation is the same as before with the simple solution
	({{T[t], =, (Subscript[T, 0]) (E^(-λ k t))}}).
Be brave and write down the other equation
	({{Δϕ[r,θ], =, -λ ϕ[r,θ]}})
This is an eigenvalue problem for the Laplacian. 
Be very brave and plug a separated solution ϕ[r,θ]=R[r]Θ[θ] into the expression for the polar Laplacian to get
	(1/(r^2))R[r]Θ''[θ]+ (1/r)(d/dr)(r R'[r])Θ[θ]=-λ R[r] Θ[θ]
Divide by R[r]  and Θ[θ] then mulitply by (r^2) to get
	((Θ''[θ])/(Θ[θ]))+ (r/( R[r]))(d/dr)(r R'[r])=-λ (r^2) 
and put the R stuff on the far side with the -λ
	((Θ''[θ])/(Θ[θ]))=-λ (r^2)-(r/(R[r]))(d/dr)(r R'[r]) 
argue (just like before) that this can only work if the two sides are equal to a separation constant (we will call it -μ because λ is already taken) and separate to get
	((Θ''[θ])/(Θ[θ]))=-μ 
and 
	-λ (r^2)-(r/(R[r]))(d/dr)(r R'[r])=-μ .
As always, we do the simplest equation first
	Θ''[θ]+μ Θ[θ]=0 
which has wiggly solutions iff μ=(ω^2)>0
	(Subscript[C, 1])Cos[ω θ]+(Subscript[C, 2])Sin[ω θ]
and it also has a constant solution if μ=0. These wiggly solutions satisfies periodic conditions if ω=n so we need
	μ=(n^2).
All we have to do now is plug stuff back in to the only remaining equation (the one for R) and sort out what it is telling us. First just plug the separation constant back in.
	-λ (r^2)-(r/( R[r]))(d/dr)(r R'[r])=-(n^2) 
then tidy up by multiplying by R
	-λ (r^2)R[r]-r(d/dr)(r R'[r])=-(n^2) R[r]
to get the ODE
	-λ (r^2)R[r]-r (r R''[r] +R'[r])=-(n^2) R[r]
or the slightly cleaner
	(r^2) R'' + r R'+(λ (r^2)-(n^2))R=0.
This is a Bessel equation and the solution is
	(Subscript[c, 1]) (BesselJ[n, r (Sqrt[λ])])+(Subscript[c, 2]) (BesselY[n, r (Sqrt[λ])])
where  (BesselJ[n, r (Sqrt[λ])])  and (BesselY[n, r (Sqrt[λ])]) are Bessel functions of the first and second kind.

R[r]/.DSolve[ (r^2)R''[r]+r R'[r]+(λ (r^2)-(n^2))R[r]==0,
R[r],
r]

Just like the “0” case we saw before the Js are smooth while (lim_(r->(SuperPlus[0])))(Subscript[Y, n])(r)=-∞.

Plot[
Evaluate[
Table[
BesselJ[n,r],
{n,0, 4}]
],
{r, 0, 8},
PlotLegends->Range[0,4]]

Plot[
Evaluate[
Table[
BesselY[n,r],
{n,0, 4}]
],
{r, 0, 6},
PlotLegends->Range[0,4],
PlotRange->{-1,0.7}]

ParametricPlot3D[
{r Cos[θ], r Sin[θ], BesselJ[0,r  BesselJZero[0,1]]},
{r,0,1},
{θ,0,2π}]

ParametricPlot3D[
{r Cos[θ], r Sin[θ], BesselJ[0,r  BesselJZero[0,2]]},
{r,0,1},
{θ,0,2π}]

ParametricPlot3D[
{r Cos[θ], r Sin[θ], Cos[θ] BesselJ[1,r  BesselJZero[1,1]]},
{r,0,1},
{θ,0,2π}]

ParametricPlot3D[
{r Cos[θ], r Sin[θ], Cos[2θ] BesselJ[2,r  BesselJZero[2,1]]},
{r,0,1},
{θ,0,2π}]

ParametricPlot3D[
{r Cos[θ], r Sin[θ], Cos[2θ] BesselJ[2,r  BesselJZero[2,2]]},
{r,0,1},
{θ,0,2π}]

kMax=4;
Plot[
Evaluate[
Table[
BesselJ[0,BesselJZero[1,k] r],
{k,1, kMax}]
], 
{r,0, 1},
PlotLabel->"Neumann"]

Big picture: words

We look for special solutions of the form u[r,t]=R[r]T[t] by plugging this ansatz into the PDE and BC and separating the variables.

We argue that the two parts of the separated equation have to equal the same constant λ.

Solving the “r” ODE for R[r] and satisfying the BC gives a discrete set of separation constants λ.

Plugging those values for λ into the “t” ODE for T[t] gives the time dependence for each of the special solutions.

Weighted sums of these special solutions are also solutions since the PDE is linear.

Turns out the weighted sums are sufficiently general that every solution can be represented in that form.

The time dependence of the weighted sums gives us really interesting information about the solutions.

(Subscript[u, t])=Δu on Ω  with u(Subscript[|, ∂ω])=0

We are going to think about and compare the unit circle and the unit interval.  The time dependence is exponential for both geometries. The special shapes are eigenfunctions for the Laplacian.

Clear[X,x,R,r,n,T,t]
X[n_][x_]:= Sin[n π x]
R[n_][r_]:=R[n][r]= BesselJ[0,r BesselJZero[0,n]]
T[t_,λ_]:=(E^(-λ t))
FullSimplify[{
((Laplacian[X[n][x],{x},"Cartesian"])/(X[n][x])),
((Laplacian[R[n][r],{r,θ},"Polar"])/(R[n][r]))
}]

Manipulate[
Plot3D[  
r=(Sqrt[(x^2)+(y^2)]);
R[n][r],
Manipulate[
Plot3D[  
r=(Sqrt[(x^2)+(y^2)]);
R[n][r],
{x,-1,1},
{y,-1,1},
RegionFunction->Function[{x,y,z},(x^2)+(y^2)<1]
],
{n,{1,2,3}}
]

],
{n,{1,2,3}}
]

Manipulate[
Plot3D[  
r=(Sqrt[(x^2)+(y^2)]);
R[n][r],
{x,y}∈Disk[{0,0},1]
],
{n,{1,2,3}}
]

Manipulate[
ParametricPlot3D[ {r Cos[θ],r Sin[θ], R[n][r]},
{r,0,1},
{θ,0, 2π}],
{n,{1,2,3}}
]

(Subscript[u, t])=k Δu with u(Subscript[|, ∂ω])=0
0<r<rMax and 0<θ<2π

We have three variables θ,r, and t.  This is our first three variable problem.  The separation trick still works it just gets a little messier. 

As before assume that u=T[t] ϕ[r,θ] and plug in to the PDE 
	({{T'[t]ϕ[r,θ], =, k T[t] Δϕ[r,θ]}})
and separate to get
	({{((T'[t])/(k T[t])), =, ((Δϕ[r,θ])/(ϕ[r,θ]))}})=-λ
The time equation is the same as before with the simple solution
	({{T[t], =, (Subscript[T, 0]) (E^(-λ k t))}}).
Be brave and write down the other equation
	({{Δϕ[r,θ], =, -λ ϕ[r,θ]}})
This is an eigenvalue problem for the Laplacian. 
Be very brave and plug a separated solution ϕ[r,θ]=R[r]Θ[θ] into the expression for the polar Laplacian to get
	(1/(r^2))R[r]Θ''[θ]+ (1/r)(d/dr)(r R'[r])Θ[θ]=-λ R[r] Θ[θ]
Divide by R[r]  and Θ[θ] then mulitply by (r^2) to get
	((Θ''[θ])/(Θ[θ]))+ (r/( R[r]))(d/dr)(r R'[r])=-λ (r^2) 
and put the R stuff on the far side with the -λ
	((Θ''[θ])/(Θ[θ]))=-λ (r^2)-(1/(r R[r]))(d/dr)(r R'[r]) 
argue (just like before) that this can only work if the two sides are equal to a separation constant (we will call it -μ because λ is already taken) and separate to get
	((Θ''[θ])/(Θ[θ]))=-μ 
and 
	-λ (r^2)-(r/(R[r]))(d/dr)(r R'[r])=-μ .
As always, we do the simplest equation first
	Θ''[θ]+μ Θ[θ]=0 
which has wiggly solutions iff μ=(ω^2)>0
	(Subscript[C, 1])Cos[ω θ]+(Subscript[C, 2])Sin[ω θ]
and it also has a constant solution if μ=0. These wiggly solutions satisfies periodic conditions if ω=n so we need
	μ=(n^2).
All we have to do now is plug stuff back in to the only remaining equation (the one for R) and sort out what it is telling us. First just plug the separation constant back in.
	-λ (r^2)-(r/( R[r]))(d/dr)(r R'[r])=-(n^2) 
then tidy up by muliplying by R
	-λ (r^2)R[r]-r(d/dr)(r R'[r])=-(n^2) R[r]
to get the ODE
	-λ (r^2)R[r]-r (r R''[r] +R'[r])=-(n^2) R[r]
or the slightly cleaner
	(r^2) R'' + r R'[r] +(λ (r^2)-(n^2))R=0.
This is a Bessel equation and the solution is
	(Subscript[c, 1]) (BesselJ[n, r (Sqrt[λ])])+(Subscript[c, 2]) (BesselY[n, r (Sqrt[λ])])
where  (BesselJ[n, r (Sqrt[λ])])  and (BesselY[n, r (Sqrt[λ])]) are Bessel functions of the first and second kind.

R[r]/.DSolve[ (r^2)R''[r]+r R'[r]+(λ (r^2)-(n^2))R[r]==0,
R[r],
r]

Just like the “0” case we saw before the Js are smooth while (lim_(r->(SuperPlus[0])))(Subscript[Y, n])(r)=-∞.

Now just like before the zeros fo these functions are known so we can work out the values of λ that make these functions zero at r=1 to enforce our boundary conditions.

TableForm[
Table[
N[BesselJZero[n,k]],
{n,0, 5},
{k,1,4}],
TableHeadings->{
Table[StringForm["n=``",n],{n,0,5}],
Table[StringForm["k=``",k],{k,1,4}]}
]

Our shape functions have two parameters “n” is the number of zeros in the θ direction and “k” is the number of humps in the r direction.

TabView[
Table[
StringForm["n=``",n,k]->TabView[
Table[
γ=N[BesselJZero[n,k]];StringForm["k=``",n,k]->ParametricPlot3D[
{r Cos[θ], r Sin[θ], Cos[n θ] BesselJ[n, r γ]},
{r,0, 1},
{θ,0,2π}],
{k,1,5}]
],
{n,0,4}]
]

Pic Maker

Plot[Evaluate[
 Table[
BesselJ[n,r],
{n,0,7}]
],
{r,0, 10},
PlotLegends->Table[
StringForm["(BesselJ[``, r]) ",n],{n,0,7}]
 ]

Plot[Evaluate[
 Table[
BesselY[n,r],
{n,0,7}]
],
{r,0, 10},
PlotRange->{-1.5,1},
PlotLegends->Table[
StringForm["(BesselK[``, r]) ",n],{n,0,7}]
 ]