function [x,y,V,charge,totalCharge]=laplaceOptimPlates

tic;
xgrid=200;
ygrid=200;

%create x,y mesh
x=linspace(-1,1,xgrid);
y=linspace(-1,1,ygrid);
[X,Y]=meshgrid(x,y);

%set boundary conditions
leftBC=0;
rightBC=0;
topBC=0;
bottomBC=0;

V=zeros(xgrid,ygrid); %creates a potential matrix of all zeros
%V=1-sqrt(X.^2+Y.^2);


%set the appropriate boundary conditions
V(1,:)=leftBC;
V(end,:)=rightBC;
V(:,1)=bottomBC;
V(:,end)=topBC;

%define the location of the left plate
xL1=round(xgrid*.35);
xL2=round(xgrid*.65);
yL1=round(ygrid*.45);
yL2=round(ygrid*.451);

%define the location of the right plate
xR1=round(xgrid*.35);
xR2=round(xgrid*.65);
yR1=round(ygrid*.55);
yR2=round(ygrid*.551);

maxIterations=2000; %The maximum number of times the iteration procedure will run

change=0;

%can use waitbar to show progress, but it increases runtime dramatically
%h=waitbar(0, sprintf('%1.0e',change),'Name', 'Change is ');

for i=1:maxIterations;
    Vold=V; %
    
    %establish V matrices shifted in each direction needed for vectorization
    Vdown=circshift(V,[1,0]);
    Vup=circshift(V,[-1,0]);
    Vright=circshift(V,[0,1]);
    Vleft=circshift(V,[0,-1]);
    
    %This averaging step solves Laplace's equation
    V=(Vdown+Vup+Vright+Vleft)/4;

    %re-establish the set potential of the plates
    V(xL1:xL2,yL1:yL2)=1; 
    V(xR1:xR2,yR1:yR2)=-1;
    
    %maintain the boundary conditions
    V(1,:)=leftBC;
    V(end,:)=rightBC;
    V(:,1)=bottomBC;
    V(:,end)=topBC;
    
    %track the change in the matrix at each step to know when to quit
    deltaV=V-Vold;
    change=abs(sum(sum(deltaV))/(xgrid*ygrid));
    
    %waitbar(i/maxIterations,h,sprintf('%1.0e',change));
    
    %conditional to quit the method when changes are within tolerance
    if i>1000 && change<1e-11
        break
    end
end

disp('Number of iterations')
disp(i);

disp('change parameter at final step')
disp(change);


%close(h)
figure
meshc(X,Y,V);


%Find the electric field as the gradient of the potential
[EX,EY]=gradient(-V,2/xgrid);

%Plot the vector electric field
figure
quiver(X,Y,EX,EY);

%Use Poisson's equation to find the charge density
charge=-del2(V,2/xgrid);

%plot the charge density
figure
meshc(X,Y,charge);

%set all charge values less than 1% of the maximum to zero
charge(abs(charge)<max(charge(:)/100))=0;

%sum the absolute value of the charges
totalCharge=sum(abs(charge(:)));

toc;

end

