% assuming you have made a matrix K which is the cotangent Laplacian
% for a mesh of the unit disk with vertices p and triangles t and saved it
load 'kdisk.mat'
%look at sparcity pattern of K
spy(K)

% The radii. Could use cart2pol.
rads = (p(:,1).^2 + p(:,2).^2).^(1/2);
% a crude way to find a vertex in the middle
[m,c]=min(rads);

%source in middle 
b=zeros(size(p,1),1);
b(c)=1;

%Solve Neumman zero boundary value problem for this source 
u = K(:,2:end)\b;
u=[0;u]; % zero arbitrarily at node 1 as one extra condition needed

%Contour plot. Not a particularly good way to do it but a start you could
% improve on!
F= TriScatteredInterp(p(:,1),p(:,2),u);
xi = -1:0.1:1;
[qx,qy]=meshgrid(xi,xi);
qz=F(qx,qy);
contour(qx,qy,qz)

figure
% Trisurf. Not a very good way to visualize it either!
trisurf(t,p(:,1),p(:,2),u)

%boundary nodes, as it is a disk they are easy to find, but this
% is a quick and dirty way to do it and not very robust. A more proper way 
% is to make a list of edges which only lie in one triangle and then list
% the vertices in thos edges
bvtx = find(rads >=0.99);


%Dirichlet  problem
% choosing a problem where we know the solution as it is 
% a haronic function.
f = p(bvtx,1).^2 - p(bvtx,2).^2;  % x^2- y^2 at boundary
% list of interior vertices
ivtx = setdiff(1:size(p,1),bvtx)';
% You should think about this and convince yourself it works
% it helps me to draw a picture of the matrices
uu = K(:,ivtx ) \(-K(:,bvtx)*f);
% so now we have the solution on the interior vertices.
% Assemble it in to a vector of teh solution at all vertices.
U(ivtx)=uu;
U(bvtx)=f;

%Plot. Not very well, but good enough to see it worked!

 F= TriScatteredInterp(p(:,1),p(:,2),U');
qz=F(qx,qy);
contour(qx,qy,qz)


