The rule for computational questions like this are that you can do absolutely whatever you want, as long as you get an answer to the question that you can reasonably justify. You may want to do more in order to understand the problem better. Remember that in the end the problem is purely for your own benefit.
If you check out Lemma 1.41-1.42 on page 23 of Shimura's "Introduction to ... " book
http://modular.fas.harvard.edu/edu/Fall2003/252/references/Shimura-Intro/index.html
you'll find a way to describe the cusps for Gamma(N). Two cusps a/b and c/d (in reduced form) are equivalent if and only if +/-[a,b] = [c,d] (mod N). Also, the number of Gamma(N) cusps is N^2/2*prod_{p|N}(1-p^(-2)) (for N>2). I think this means the cusps for Gamma(N) are in bijection with the following set:
Take vector +/- [a,b] where a and b are numbers mod N and GCD(a,b,N) = 1, and the +/- means that we identify [a,b] with [-a,-b]. Then given such a vector, there exists integers a',b' with gcd(a',b')=1, and a'=a(mod N) and b'=b(mod N), and the corresponding cusp is a'/b'.
Once you have the cusps for Gamma(N), you can get the cusps for Gamma1(N) and Gamma0(N) by taking orbits for the actions of the corresponding subgroups of SL_2(Z/NZ).
I have no idea if this is what MAGMA does. Algorithmically, if you have a test for cusp equivalence (and there are simple such tests), you can just choose small random rational numbers and usually generate up all cusps.