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]]
```

I do not have sufficient knowledge, to answer this question, I will delete my answer . as I deleted my answer (it has also deleted your comment) you'll have to re-ask your question about:

Sorry Arnab.