function [HarVec]= Harmonic(Matrix, i, j)
%The input Matrix is the Laplacian matrix of some graph. The scalar inputs
%i and j correspond to fixed vertices in the matrix. The function attempts to
%find the harmonic function associated with these vertices. The harmonic
%function, HarVec, must satisfy:
%1. HarVec(i)=1
%2. HarVec(j)=0
%3. Matrix*HarVec=I, where I is a vector with all zero entries except for
% entry corresponding to the vertex i.
[m,n]=size(Matrix);
%Since the i'th and j'th values of HarVec are unknown, we can get rid of
%the i'th and j'th linear equations coming from condition 3. In
%particular, the i'th and j'th rows are set to zero, except for the
%diagonal entries, which are set to one.
Matrix(i,:)= 0;
Matrix(j,:)=0;
Matrix(i,i)=1;
Matrix(j,j)=1;
I=zeros(m,1);
I(i)=1;
%I is a vector as defined earlier. The i'th entry is now one, because of
%our earlier changes.
%Now we can use matrix operations to solve for HarVec.
C=Matrix^(-1);
HarVec=C*I;
end