x, y = var('x y')
e1 = x^2 + y^2 == 41
e2 = y == x + 1
show(e1)
show(e2)

s = solve([e1, e2], [x, y])
show(s)

s1 = solve(e1, y); show(s1)

y1 = s1[0].rhs(); show(y1)

y2 = s1[1].rhs(); show(y2)

set_verbose(-1)
plot([y1, y2, e2.rhs()], (x, -10, 10), aspect_ratio = 1)

