(linear form in x and y) >= 0
the subset S is what is called a "polyhedron" in R^2.
The problem in your original post can now be
rephrased as:
    Find all integral points in the polyhedron S.
An introduction to polyhedra in Sage is at:
The polyhedron S can be input as
    S = Polyhedron(ieqs=[[-1, 1, 0], [9, -1, 0], [-1, 0, 1], [9, 0, -1]], eqns=[[-15, 1, 1]]), 
Check that our input represents the correct polyhedron:
    sage: print(S.Hrepresentation_str())
    x0 + x1 ==  15
        -x0 >= -9
         x0 >=  6
Find all integral points:
    sage: S.integral_points()
    ((6, 9), (7, 8), (8, 7), (9, 6))