GB for I=ideal(x^3-3*x^2-y+1,-x^2+y^2-1) leaves out a further step.

savecancel

Whenever the monomial/term ordering is fixed, the reduced Groebner basis of an ideal is unique (up to reordering of the list of elements). In GAP you used the lexicographic ordering, and in SageMath you didn't specify a monomial ordering, so it used the default of degrevlex, which is not lex, so you got a different result. Here is how to specify the monomial ordering (and get the same result as in GAP, up to reordering):

sage: R.<x, y> = PolynomialRing(QQ, order='lex')
sage: I = R.ideal(x^3 - 3*x^2 - y + 1, -x^2 + y^2 - 1)
sage: I.groebner_basis()
[x^2 - y^2 + 1, x*y - x - y^4 + 11*y^2 + 3*y - 13, y^5 + y^4 - 11*y^3 - 17*y^2 + 9*y + 17]

Or conversely, use GAP to get the same result as in SageMath, by choosing the degrevlex (a.k.a. grevlex) monomial ordering:

gap> R := PolynomialRing(Rationals, [ "x", "y" ]);;
gap> x := IndeterminatesOfPolynomialRing(R)[1];;
gap> y := IndeterminatesOfPolynomialRing(R)[2];;
gap> degrevlex := MonomialGrevlexOrdering(x,y);;
gap> ideal := [x^3-3*x^2-y+1, -x^2+y^2-1];;
gap> G := ReducedGroebnerBasis(ideal, degrevlex);;
gap> Display(G);
[ x^2-y^2+1, x*y^2-3*y^2-x-y+4, y^4-x*y-11*y^2+x-3*y+13 ]