You can do this in Sage but "in disguise". The idea is that if $f(z)$ is a cusp form for $\Gamma(N)$, then $g(z) := f(Nz)$ is a cusp form for a certain subgroup intermediate between $\Gamma_0(N^2)$ and $\Gamma_1(N^2)$ which Sage calls $\Gamma_H(N^2, [N + 1])$; the $N + 1$ is here because it generates the subgroup of $(\mathbf{Z} / N^2 \mathbf{Z})^\times$ consisting of classes that are 1 mod $N$. If $f(z) = \sum a_n q^{n/N}$, then $g(z) = \sum a_n q^n$, so if you have the $q$-expansion of $g$ then you have the $q$-expansion of $f$ and vice versa.

```
----------------------------------------------------------------------
| Sage Version 5.9, Release Date: 2013-04-30 |
| Type "notebook()" for the browser-based notebook interface. |
| Type "help()" for help. |
----------------------------------------------------------------------
sage: G = GammaH(36, [7])
sage: G.index() == Gamma(6).index()
True
sage: [g.q_expansion(25) for g in CuspForms(G, 2).basis()]
[q - 4*q^7 + 2*q^13 + 8*q^19 + O(q^25)]
```

(Note that this $g$ actually has trivial character, and is of CM type, so its $q$-expansion has lots of zero coefficients.)