%% Plot of the solution using Fourier modes
% This code computes the solution of the problem 2.5.9 (b) in your 
% texboook. Note that only a partial sum of the actual series solution can be plotted.  

a=0.4;  b=1; % bounds for r
nr=50; nt=50; % discretization of the domain (in poloar coodrinates)

[R,T] = meshgrid(a:(b-a)/nr:b,0:pi/(2*nt):pi/2);

N=20; % number of Fourier modes

w=zeros(N,1);
A=zeros(N,1);
U=zeros(nr+1,nt+1);

L=log(b/a);

for n=1:N
    w(n)=n*pi/L;
    
    f=@(r) (1-r.^2).*sin(w(n).*log(r/a))./r;
    A(n)=1/sinh(w(n)*pi/2)*2/L*quad(f, a, b);
    
    U=U+A(n)*sin(w(n)*log(R/a)).*sinh(w(n)*T);
end

X=R.*cos(T); Y=R.*sin(T); % change of variables from polar to rectangular

subplot(1,2,1),  surf(X,Y,U), axis square
subplot(1,2,2), surf(X,Y,U), axis square, view(2)