function evalGH2 (loc,k)
% Recursive function to evaluate g(k,k) and h(k,k) for current busy period
% Uses backtracking to expand g() and h() for ancestors as necessary
global cache
% By assumption, k-arrival cache entries do not exist at the current node
v= k-cache(loc).i;
if k==1 % Alg 1: lines 9,10
    cache(loc).g(atZero(0))=cache(loc).size;
    cache(loc).h(atZero(0))=cache(loc).size^2/2;
else if cache(loc).i ==1 % Alg 1: lines 12, 13
        cache(loc).g(atZero(v))=cache(loc).g(atZero(v-1))*cache(loc).size/k;
        cache(loc).h(atZero(v))=cache(loc).g(atZero(v-1))*cache(loc).size^2/2;
    else % Attempt to execute last iteration of for loop, lines 14 - 24
        up=cache(loc).parent;
        if numel(cache(up).g) <atZero(v+1)
            evalGH2(up,k); % backtrack to fill cache from previous iteration
        end
        cache(loc).g(atZero(v)) = 0; % line 15
        if cache(loc).i>k-1 % min function on line 16
            cache(loc).h(atZero(v))=cache(up).g(atZero(v));
        else
            cache(loc).h(atZero(v))=cache(loc).g(atZero(v-1));
        end
        cache(loc).h(atZero(v)) = cache(loc).h(atZero(v)) * ...
            cache(loc).size*(cache(up).X + cache(loc).size/2);
        top=v+1;
  %      f=1;
        cache(loc).f(atZero(top)) = cache(loc).f(atZero(top-1)) * cache(loc).size/top;
        for m=0:top % Convolution in lines 20 - 23
            cache(loc).g(atZero(v))=cache(loc).g(atZero(v))+...
                cache(loc).f(atZero(m))*cache(up).g(atZero(v-m+1));
            cache(loc).h(atZero(v))=cache(loc).h(atZero(v))+...
                cache(loc).f(atZero(m))*cache(up).h(atZero(v-m+1));
        end
    end
end