We present a new numerical method for computing the GSVD of two matrices A and B. This method is a variation on Paige's method. It differs from previous algorithms in guaranteeing both backward stability and convergence. There are two innovations. The first is a new preprocessing step which reduces A and B to upper triangular forms satisfying certain rank conditions. The second is a new 2 by 2 triangular GSVD algorithm, which constitutes the inner loop of Paige's method. We present proofs of stability and convergence of our method, and demonstrate examples on which all previous algorithms fail.