The statement of the problem (Subsection 2.5.2 in the book):

In this file I will consider the Laplace's equation in a disk. See Subsection 2.5.2 (page 73) in the book. The equation is

(∇^2)u=(1/r)(∂/(∂r)) (r ((∂u)/(∂r)))+(1/(r^2))(((∂^2)u)/(∂(θ^2)))=0,   0 < r≤R,    -π ≤ θ ≤π,   
  θ   "boundary conditions"   {({{u(r,-π)-u(r,π)= 0,}, {((∂u)/(∂θ))(r,-π)-((∂u)/(∂θ))u(r,π)= 0,}})
  r     boundary conditions    {({{u(R,θ)=f(θ),    0 ≤θ≤2π,}, {|u(0,θ)| < ∞ . }})

We will solve this boundary value problem by the separation of variables method. We look for the solution of the form   u(r,θ)=  A(r) B(θ).  This leads to :

(∇^2)u=((1/r)(d/(d r)) (r ((d )/(d r))A(r))) B(θ)+(1/(r^2))A(r)((d^2)/(d (θ^2)))B(θ)=0,
    B(-π) - B(π) = 0,
B'(-π) - B'(π) = 0 .

Here, as before, we ignore the nonhomogeneous set of boundary conditions.  Separating variables we obtain 

((r(d/(d r)) (r ((d )/(d r))A(r)) )/(A(r)))=-((((d^2)/(d (θ^2)))B(θ))/(B(θ)))=λ, 

what leads to the boundary eigenvalue problem for the function B(θ) :

-(((d )^(2 ))/(d (θ^2)))B(θ) = λ B(θ),  
       B(-π) - B(π)  = 0,         
    B'(-π) - B'(π) = 0 ,

and the equation

r(d/(d r)) (r ((d )/(d r))A(r)) - λ A(r)= 0,

for the function A.  

The boundary eigenvalue problem for the function A is identical to the problem studied in 2.5.1 (thin circular ring).  The eigenvalues of this problem and the corresponding eigenfunctions are given by: 

The eigenvalues:      The corresponding eigenfunction(s): 

  (Subscript[λ, 0]) =  0 ,                 1  (the constant function), 
  (Subscript[λ, 1]) =  (1^(2 )),               Sin[x],   and  Cos[x]  (two linearly independent eigenfunctions),  
  (Subscript[λ, 2]) =  (2^(2 )),               Sin[2x],   and  Cos[2x]  (two linearly independent eigenfunctions),   
  
  and in general: 
  
  (Subscript[λ, n]) =  (n^(2 )),               Sin[n x],   and  Cos[n x] , n=0,1,2, ...  . 
  
Note that the last formula makes sense for n=0, since Cos[0 x] = 1, and Sin[0 x] = 0 (which can be ignored).

Solving the equation for the function A:

We conclude that the solutions for the boundary eigenvalue problem for the function B are numbers λ=  (n^(2 )) (eigenvalues) and the corresponding functions Cos[n x],Sin[n x], n=0, 1,2, 3, ... (eigenfunctions).  
For these values λ = (n^(2 )) we have to solve the equation for the function A:   

 
     r(d/(d r)) (r ((d )/(d r))A(r)) - (n^2) A(r)= 0 ,    
  

or in a more transparent form    

 
     (r^2) ((d^(2 ))/(d (r^2)))A +r ((d )/(d r))A- (n^2) A= 0  .   



Let us show how to use Mathematica to solve this equation. The basic command is

DSolve[(r^2)A''[r]+r A'[r]-(n^2)A[r]==0,A[r],r]

But this is not a desirable form of the solution; there must be a simpler expression for

TrigToExp[{Cosh[n Log[r]],Sinh[n Log[r]]}]

Since any linear combinations of solutions is a solution we have that

{{1,1},{1,-1}}.{((r^(-n))/2)+((r^n)/2),-((r^(-n))/2)+((r^n)/2)}

are also solutions. Thus, Mathematica should have given the general solution

C[1](r^n) + C[2](r^(-n))

Below, just to demonstrate some Mathematica commands, I will try to force Mathematica to give this solution.

Collect[Expand[Simplify[TrigToExp[DSolve[(r^2)A''[r]+r A'[r]-(n^2)A[r]==0,A[r],r][[1]][[1]][[2]]],r>0]],(r^n)]

In the last expression constants are complex numbers and we can replace them with arbitrary constants C[1] and C[2] and we can write this general solution as:

solA=Collect[Expand[Simplify[TrigToExp[DSolve[(r^2)A''[r]+r A'[r]-(n^2)A[r]==0,A[r],r][[1]][[1]][[2]]],r>0]],(r^n)]/.{(-(1/2) I C[2]+((C[1])/2))-> C[1],((1/2) I C[2]+((C[1])/2))-> C[2]}

Note that the above solution is the general solution of the equation for A only if n > 0. Namely, for the above formula to be the general solution the functions (r^n) and  (r^(-n)) must be linearly independent on 0 ⩽r⩽ R.  For n = 0 this obviously is not a case.  Thus we have to solve the equation for n = 0  separately:

DSolve[(r^2)A''[r]+r A'[r]==0,A[r],r][[1]][[1]][[2]]

Nice solution!  Now we recall the boundary condition for u(r,θ) near r=0 : 

|u(0,θ)| < ∞ , or in terms of A :  |A(near 0)| < ∞. 

This clearly eliminates the functions  (r^(-n))  and  Log[r] as possible solutions for A(r).  Thus the only possible solutions for A are 

(Subscript[A, n])(r) =  C[1] (r^n) ,  n = 0,1,2, ... , 

where C[1] is and arbitrary constant.

Conclusion (the solution of Laplace's equation):

We have found infinitely many solutions for the Laplace's equation on a circular disk: 

((r^n)/(R^n))Sin[n θ] ,     ((r^n)/(R^n))Cos[n θ]   where n=0,1,2,3, ... .  

Using the principle of superposition we conclude that the function

   u(r,θ)=(Subscript[a, 0])+(∑_(n=1)^∞)(Subscript[a, n])((r^n)/(R^n))Cos[n θ]  + (∑_(n=1)^∞)(Subscript[b, n])((r^n)/(R^n))Sin[n θ]      

is also a solution.  We just need to determine the constants (Subscript[a, n]) and (Subscript[b, n]) in such a way that the boundary conditions 

 u(R,θ) = f(θ), 0 ≤ θ ≤2π  

is satisfied.  Put r=R in the above formula for u(r,θ) and we get 


  
   f(θ)= u(R,θ) =(Subscript[a, 0])+(∑_(n=1)^∞)(Subscript[a, n])Cos[n θ]  + (∑_(n=1)^∞)(Subscript[b, n])Sin[n θ]   
				. 

The orthogonality of  the functions Sin[n x] and Cos[n x] leads to the formula for (Subscript[a, m]):

  (∫_(-π)^π)f(θ) Cos[m θ] dx
		=(Subscript[a, m])(∫_(-π)^π)Cos[m θ] *Cos[m θ] dθ. 

Thus:
  (Subscript[a, m])= ((  (∫_(-π)^π)f(θ) Cos[m θ] dθ)/((∫_(-π)^π)Cos[m θ] *Cos[m θ] dθ)), m =0,1,2,3, ... . 

Similarly: 

  (Subscript[b, m])= ((  (∫_(-π)^π)f(θ) Sin[m θ] dθ)/((∫_(-π)^π)Sin[m θ] *Sin[m θ] dθ)), m =1,2,3, ... . 

Since we can calculate:

FullSimplify[{(∫_(-π)^π)Cos[0 x] *Cos[0 x] \[DifferentialD]x,(∫_(-π)^π)Cos[n x] *Cos[n x] \[DifferentialD]x,

(∫_(-π)^π)Sin[n x] *Sin[n x] \[DifferentialD]x},n∈Integers]

we can rewrite our formulas for (Subscript[a, n]) and (Subscript[b, n]) as 

  (Subscript[a, 0])= ((  1)/(2 π))(∫_(-π)^π)f(θ) dθ  and  

  (Subscript[a, n])= ((  1)/π) (∫_(-π)^π)f(θ) Cos[n θ] dθ, n =1,2,3, ...  and  

  (Subscript[b, n])= ((  1)/π) (∫_(-π)^π)f(θ) Sin[n x] dθ, n =1,2,3, ... . 

Thus the solution of the heat equation for a circular disk is: 

  u(r,θ)=(Subscript[a, 0])+(∑_(n=1)^∞)(Subscript[a, n])((r^n)/(R^n))Cos[n θ]  + (∑_(n=1)^∞)(Subscript[b, n])((r^n)/(R^n))Sin[n θ]   

with (Subscript[a, n]) and (Subscript[b, n]) given by the above formulas.

Next we implement these formulas in Mathematica:

Mathematica implementation of the solution:

Here is our function f(θ)  which we will approximate with the first 15 (or nn) terms.

Clear[ff];
ff[θ_]=((((θ^2)-(π^2)))^2)(θ+π-2)Exp[-θ-π]

test={(ff[θ]/.{θ->-π}),(ff[θ]/.{θ-> π}),(D[ff[θ],θ]/.{θ-> -π}),(D[ff[θ],θ]/.{θ-> π})}

Plot[{ff[θ]},{θ,-π,π}]

Clear[rR,nn,las,lbs];

rR=1;
nn=15;

las=Table[(1/π)*NIntegrate[Expand[ff[θ]*Cos[n θ]],{θ,-π,π},MaxRecursion->20,
PrecisionGoal->12,WorkingPrecision->20,AccuracyGoal-> 12],
		{n,1,nn}];
lbs=Table[(1/π)*NIntegrate[Expand[ff[θ]*Sin[n θ]],{θ,-π,π},MaxRecursion->20,
PrecisionGoal->12,WorkingPrecision->20,AccuracyGoal-> 12],
		{n,1,nn}];

las

lbs

Clear[uu];

uu[r_,θ_]=(1/(2π ))NIntegrate[ff[θ],{θ,-π,π},MaxRecursion->20,
PrecisionGoal->12,WorkingPrecision->20,AccuracyGoal-> 12]+
		Sum[las[[n]]*((r^n)/(rR^n))*Cos[n θ],{n,1,nn}]+
		Sum[lbs[[n]]*((r^n)/(rR^n))*Sin[n θ],{n,1,nn}] ;

uu[.5,π/2]

uu[0,0]

Test of the approximation:

Plot[{ff[θ],uu[rR,θ]},{θ,-π,π},PlotStyle->{{Thickness[0.007],Blue},{Thickness[0.004],Red}}]

I find  maximum and minimum of ff to set the range of plots correctly.

FindMaximum[ff[θ],{θ,-1}]

FindMinimum[ff[θ],{θ,-3}]

I experimented and found out that the following view point works well for the plot below.

vp={-1.9533572861214523`,-2.405091682993332`,1.3602486203457567`}

pluu=ParametricPlot3D[ {r Cos[θ],r Sin[θ],uu[r,θ]},
               {r, 0,rR}, {θ, -π, π},Mesh->False,PlotRange->{{-1,1},{-1,1},{-11,5}},BoxRatios-> {1,1,1},PlotPoints->{20,50},ImageSize->300,ViewPoint-> Dynamic[vp]]

I want to show the unit circle on the graph as well. So, I define it separately. Below is the unit circle placed at the level below the minimum temperature.

uc=ParametricPlot3D[{rR Cos[θ],rR Sin[θ],-10.9},  {θ, -Pi, 2Pi},PlotStyle-> Thickness[0.005]]

ic=ParametricPlot3D[{rR Cos[θ],rR Sin[θ],ff[θ]},  {θ, -Pi, Pi},PlotStyle-> {Thickness[0.01],RGBColor[0,1,0]},PlotPoints-> 50,BoxRatios-> {1,1,1},ViewPoint-> vp]

Show[pluu,ic,uc,BoxRatios-> {1,1,1},ViewPoint-> vp]

I do not know how to do a density plot in polar coordinates in Mathematica. So, I fool Mathematica to think that it works in rectangular coordinates.  Interestingly, Mathematica has a function ArcTan as a function of two variables which will give exactly the angle θ. The only problem is that this function is not defined at {0,0}, so I define it to be 0 at the origin.

arctan[0.,0.]=0;arctan[x_,y_]:=ArcTan[x,y];DensityPlot[ uu[(Sqrt[(x^2)+(y^2)]),arctan[x, y]],
               {x,-1,1}, {y, -1,1},PlotRange->{-10.6,4.9},PlotPoints->{20,20},RegionFunction->Function[{x,y,z},Norm[{x,y}]<1.],Axes->False,Frame->True,LightingAngle-> {0,Pi/6},ImageSize->300]

DensityPlot[ uu[(Sqrt[(x^2)+(y^2)]),arctan[x,y]],
               {x,-1,1}, {y, -1,1},PlotRange->{-10.6,4.9},PlotPoints->{30,30},RegionFunction->Function[{x,y,z},Norm[{x,y}]<1.],Axes->False,Frame->True,ColorFunction-> (RGBColor[1,1-#^2,1-#^2]&),LightingAngle-> {0,Pi/3},Epilog->{Circle[]},ImageSize->300]