Define the affine variety $X = V (y − x^2, y − x + 1)$.

edit

The first example you cite is about points on a curve. Your variety is the intersection of two curves ($y = x^2$ and $y = x - 1$), which is a finite set of points (and not a curve), so you cannot use the same technique.

I assume by "affine variety" you mean affine variety defined over $\mathbb{Q}$, and by "rational points" you mean points defined over $\mathbb{Q}$. In that case, what you did later is correct:

sage: R.<x,y> = PolynomialRing(QQ)
sage: I = R.ideal(y-x^2,y-x+1)
sage: I.variety()
[]

This tells you that the set of points in the variety which are defined over $\mathbb{Q}$ (i.e. which have coordinates in $\mathbb{Q}$) is empty. You can also try to find (numerically) points defined over $\mathbb{R}$:

sage: I.variety(RR)
[]

Again it is an empty set, because the two curves don't intersect when drawn in $\mathbb{R}^2$. You can try to find (numerically) points defined over $\mathbb{C}$:

sage: I.variety(CC)
[{y: -0.500000000000000 - 0.866025403784439*I, x: 0.500000000000000 - 0.866025403784439*I},
 {y: -0.500000000000000 + 0.866025403784439*I, x: 0.500000000000000 + 0.866025403784439*I}]

Here we find two points in $\mathbb{C}^2$. If you are looking for exact answers, you can also try to find points defined over $\bar{\mathbb{Q}}$ (the algebraic closure of $\mathbb{Q}$):

sage: I.variety(QQbar)
[{y: -0.50000000000000000? - 0.866025403784439?*I, x: 0.50000000000000000? - 0.866025403784439?*I},
 {y: -0.50000000000000000? + 0.866025403784439?*I, x: 0.50000000000000000? + 0.866025403784439?*I}]

The result looks numerical, but Sage knows exactly what these points are, and you can do exact arithmetic with them. In this case you can also express them as radicals:

sage: [(P[x].radical_expression(), P[y].radical_expression()) for P in I.variety(QQbar)]
[(-1/2*I*sqrt(3) + 1/2, -1/2*I*sqrt(3) - 1/2),
 (1/2*I*sqrt(3) + 1/2, 1/2*I*sqrt(3) - 1/2)]

This is similar to what you would get if you asked Sage to solve the equations symbolically:

sage: solve(map(SR, I.gens()), map(SR, R.gens()))
[[x == -1/2*I*sqrt(3) + 1/2, y == -1/2*I*sqrt(3) - 1/2], [x == 1/2*I*sqrt(3) + 1/2, y == 1/2*I*sqrt(3) - 1/2]]