How can one find a basis for modular forms with rational Fourier coefficients?

edit

Let us work explicitly in a particular example. Let $p$ be the prime $211$. I suppose the question wants to have a basis for the space:

sage: S = CuspForms( Gamma0(p), 2 )
sage: S
Cuspidal subspace of dimension 17 
of Modular Forms space of dimension 18
for Congruence Subgroup Gamma0(211) of weight 2 over Rational Field

(Result was manually adjusted to fit in page.) Then a basis of this space is explicitly:

sage: for s in S.basis():
....:     print s.qexp(22)
....:     
q + 243/230*q^18 - 162/115*q^19 - 188/115*q^20 + 26/115*q^21 + O(q^22)
q^2 - 117/230*q^18 - 493/460*q^19 - 33/115*q^20 - 131/460*q^21 + O(q^22)
q^3 - 38/23*q^18 - 3/23*q^19 + 159/92*q^20 + 9/23*q^21 + O(q^22)
q^4 - 57/115*q^18 - 423/230*q^19 - 41/115*q^20 - 111/230*q^21 + O(q^22)
q^5 + 219/230*q^18 + 84/115*q^19 - 133/230*q^20 - 22/115*q^21 + O(q^22)
q^6 + 17/230*q^18 - 88/115*q^19 - 122/115*q^20 + 34/115*q^21 + O(q^22)
q^7 - 20/23*q^18 - 27/23*q^19 - 87/92*q^20 - 11/23*q^21 + O(q^22)
q^8 + 189/115*q^18 - 22/115*q^19 - 203/115*q^20 + 66/115*q^21 + O(q^22)
q^9 - 159/230*q^18 - 9/115*q^19 + 273/230*q^20 + 27/115*q^21 + O(q^22)
q^10 - 123/46*q^18 + 75/92*q^19 + 52/23*q^20 + 5/92*q^21 + O(q^22)
q^11 - 21/230*q^18 - 87/230*q^19 - 59/115*q^20 + 31/230*q^21 + O(q^22)
q^12 - 197/230*q^18 - 663/460*q^19 - 467/460*q^20 - 81/460*q^21 + O(q^22)
q^13 + 59/230*q^18 - 117/230*q^19 - 119/115*q^20 - 109/230*q^21 + O(q^22)
q^14 - 11/10*q^18 + 4/5*q^19 + 6/5*q^20 + 3/5*q^21 + O(q^22)
q^15 - 509/230*q^18 + 169/460*q^19 + 561/460*q^20 + 183/460*q^21 + O(q^22)
q^16 - 99/46*q^18 - 26/23*q^19 + 45/46*q^20 + 9/23*q^21 + O(q^22)
q^17 - 32/23*q^18 + 12/23*q^19 + 169/92*q^20 + 10/23*q^21 + O(q^22)

I hope this helps.

Note: Then the above space of newforms has the elements, shown with lower precision to fit in the web page:

sage: p = 211
sage: M = Newforms( Gamma0(p), 2, names='a' )
sage: for f in M:
....:     print f.qexp(4)
....:     
q + a0*q^2 + (a0 + 1)*q^3 + O(q^4)
q + a1*q^2 + (-a1^2 - a1 + 1)*q^3 + O(q^4)
q + a2*q^2 + (-a2 - 1)*q^3 + O(q^4)
q + a3*q^2 + (9/58*a3^8 + 15/58*a3^7 - 2*a3^6 - 157/58*a3^5 + 235/29*a3^4 + 222/29*a3^3 - 637/58*a3^2 - 161/29*a3 + 62/29)*q^3 + O(q^4)

So let us put the hands on the coefficients in degree $q^2$.

sage: [ list(f.qexp())[2] for f in M ]
[a0, a1, a2, a3]
sage: [ list(f.qexp())[2].parent() for f in M ]
[Number Field in a0 with defining polynomial x^2 - x - 1,
 Number Field in a1 with defining polynomial x^3 + 2*x^2 - x - 1,
 Number Field in a2 with defining polynomial x^3 - 4*x + 1,
 Number Field in a3 with defining polynomial x^9 + x^8 - 14*x^7 - 11*x^6 + 66*x^5 + 36*x^4 - 123*x^3 - 38*x^2 + 72*x + 8]

This is a complicated situation.