function [ParentList DiscoveryList] = DepthFirstSearch(EdgeListG, startV, n)

%We initialize the tree growing edge lists:
EdgeListT = []; %initially T consists of just the starting vertex startV (no edges)
%The boundary edges are all edges adjacent to startV
EdgeListGEndpts = EdgeListG(:,1:2); %indices have been removed (i.e., the third column was truncated)
[a b] = find(EdgeListGEndpts == startV); 
BdyEdges = EdgeListG(a,:);

%The boundary edges are the only edges that are initially "used"
UnusedEdges = setdiff(EdgeListG, BdyEdges, 'rows');

DiscoveryList = [startV];
ParentList = zeros(1,n);
%a zero entry in ParentList will mean the vertex has no parent "NULL"
%In the initial tree T = [startV], no vertex has a parent (and the root
%startV never will)

%Before entering into the while loop, we do the first iteration separately
%because initially EdgeListT is empty (so does not have startV)
%First Edge Selection Step (DFS)
%Need to a boundary edge with tree endpoint as LARGE as possible discovery
%number
%We initialize a vector of discovery numbers for the vertices so that
%startV has discovery number zero, and all other vertices have discovery
%number = -1
DiscoveryNumVec = -ones(1,n);
DiscoveryNumVec(startV) = 0;
%since all bdy edges have the same tree endpoint, any will do:
e = BdyEdges(1,:); %third entry is the index of edge w.r.t. EdgeListG
if e(1) == startV
    nonTreeEndpt = e(2);
else
   nonTreeEndpt = e(1);
end
DiscoveryList = [DiscoveryList nonTreeEndpt];
ParentList(nonTreeEndpt) = startV;
DiscoveryNumVec(nonTreeEndpt) = 1;
%%%Pasted from BdyEdgesUpdateInd program
%We need to delete any edge of BdyEdges that has nonTreeEndpt
%as one of its endpoints.
BdyEdgesEndpts = BdyEdges(:,1:2);
[a b] = find(BdyEdgesEndpts == nonTreeEndpt);
BdyEdges(a,:) = []; BdyEdgesEndpts(a,:) = [];
%Next, we find the unused edges that are incident to nonTreeEndpt.  The
%vector a below gives their row indices in UnusedEdges, and the vector b
%gives the columnn index(1 or 2) in which nonTreeEndpt occurs.
UnusedEdgesEndpts = UnusedEdges(:,1:2);
[a b] = find(nonTreeEndpt==UnusedEdgesEndpts);
%We add these edges to BdyEdges
BdyEdges = [BdyEdges; UnusedEdges(a,:)];
%and delete them from UnusedEdges
UnusedEdges(a,:)=[];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
EdgeListT = [EdgeListT; e];
EdgeListTEndpts= EdgeListT(:,1:2);
DiscNum = 1;

while ~isempty(BdyEdges)
    DiscNum = DiscNum + 1;
    %Next Edge Selection Step (DFS)
    %Need to find a  boundary edge whose tree endpoint discovery number is
    %as large as possible.
    EdgeListTEndpts= EdgeListT(:,1:2);
    BdyEdgesEndpts = BdyEdges(:,1:2);
    [Disc Disc_ind] = find(DiscoveryNumVec >=0); %ind = vertices that have been discovered
    BdyTreeVertices = intersect(Disc_ind, BdyEdgesEndpts(:));
    maxDiscNumBdyTreeVertex = max(DiscoveryNumVec(BdyTreeVertices));
    [c d] = find(DiscoveryNumVec(BdyTreeVertices) == maxDiscNumBdyTreeVertex);
    nextParentVertex = BdyTreeVertices(d(1));
    [a b] = find(nextParentVertex == BdyEdgesEndpts); %a = indices of bdy edges incident to maxBdyTreeVertices(1)  
    e = BdyEdges(a(1),:); %third entry is the index of edge w.r.t. EdgeListG
    if ismember(e(1),EdgeListT(:,1:2))
        nonTreeEndpt = e(2); j = 1;
    else
        nonTreeEndpt = e(1);  j = 2;
    end
    DiscoveryList = [DiscoveryList nonTreeEndpt];
    ParentList(nonTreeEndpt) = e(j);
    DiscoveryNumVec(nonTreeEndpt) = DiscNum;
    [BdyEdges UnusedEdges]= BdyEdgesUpdateInd(EdgeListT, BdyEdges, e(1:2), UnusedEdges);
    EdgeListT = [EdgeListT; e];
end
%Finally we format the outputs
ParentList2 = zeros(2,n);
ParentList2(1,:) = 1:n;
ParentList2(2,:) = ParentList;
ParentList = ParentList2';