%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%   Laplace's Equation on an annulus Example
%
%   In this code we plot a contour plot and the surface for a solution to
%   the following boundary value problem (in polar coordinates)
% 
%   \Delta u=0
%   u(x,y)=y^2, x^2+y^2=r1^2
%   u(x,y)=0,   x^2+y^2=r0^2
%   
%   Inputs:
%   1. n=number of discretization points
%   2. r1=outer radius of annulus
%   3. r0=inner radius of annulus
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Inputs
n=100;
r1=2;
r0=1;

%% Construction of Mesh
r=linspace(r0,r1,n);
th=linspace(0,2*pi,n);

[R,TH]=meshgrid(r,th);

%% Solution
u=4*(R-R.^(-1)).*sin(TH)-16/63*(R.^3-R.^(-3)).*sin(3*TH);

%% Conversion of mesh to Cartesian Coordinates
X=R.*cos(TH);
Y=R.*sin(TH);

%% Plotting Contour Plot
figure('units', 'normalized', 'outerposition', [0,0,1,1]);
contourf(X,Y,u);
axis square;
ylabel('$y$','Interpreter','latex','Fontsize',36);
xlabel('$x$','Interpreter','latex','Fontsize',36);
ax=gca;
ax.FontSize = 36;

%% Plotting Surface
figure('units', 'normalized', 'outerposition', [0,0,1,1]);
surf(X,Y,u,'edgecolor','none');
axis square;
ylabel('$y$','Interpreter','latex','Fontsize',36);
xlabel('$x$','Interpreter','latex','Fontsize',36);
zlabel('$u(x,y)$','Interpreter','latex','Fontsize',36);
ax=gca;
ax.FontSize = 36;