Computing the q-series of the j-invariant It is a fundamental fact, often quoted these days in its connection with Monstrous Moonshine, that the q-expansion (i.e., the Laurent expansion in a neighborhood of $\tau = i\infty$) of the j-invariant is given by
$$j(\tau) = \frac{1}{q} + 744 + 196884q + 21493760q^2+\dots, \quad q = e^{2 \pi i \tau}.$$
I do not recall, however, ever seeing a modern treatment, nor even a hint, of how one might go about obtaining this expansion. Does anybody know a nice way to compute these coefficients? (I mean a way which does not invoke Moonshine, not that I'd expect that to make the computation more pleasant.) Is there a standard way to do it?
I did find one approach published by H.S. Zuckerman in the late 1930s*, which makes use of a "fifth order multiplicator equation" for $j(\tau)$ -- distilled from Fricke and Klein's Vorlesungen uber die Theorie der elliptischen Modulfunktionen --  and an identity of Ramanujan for the generating function of partition numbers of the form $p(25n + 24)$. Is this typical?
*Zuckerman, Herbert S., The computation of the smaller coefficients of $J(\tau)$. Bull. Amer. Math. Soc. 45, (1939). 917–919.}
 A: The trick is to write the $j$-invariant function in terms of Eisenstein series, whose $q$-expansions have a simple expression. See Silverman's "Advanced Topics in the Arithmetic of Elliptic Curves", Chapter 1, Section 7, and in particular, Proposition 7.4 and Remark 7.4.1.
In particular 
$$
j(\tau) = 1728\frac{g_2(\tau)^3}{\Delta(\tau)}=1728 \frac{g_2(\tau)^3}{g_2(\tau)^3-27g_3(\tau)^2} = \frac{\left(1+240\sum_{n\geq 1}\sigma_3(n)q^n\right)^3}{q\prod_{n\geq 1}(1-q^n)^{24}},
$$
or, if you prefer,
$$
j(\tau) = 1728 \frac{\left(1+240\sum_{n\geq 1}\sigma_3(n)q^n\right)^3}{\left(1+240\sum_{n\geq 1} \sigma_3(n)q^n\right)^3 - \left(1-504\sum_{n\geq 1}\sigma_5(n)q^n \right)^2}.
$$
A: This is addressed in various texts, including the Silverman book cited by A.Lozano-Robledo and (if memory serves) Chapter VII of Serre's A Course in Arithmetic.  Explicitly, let $E_4$ and $E_6$ be the normalized Eisenstein series
$$
E_4 = 1 + 240 \sum_{n=1}^\infty \frac{n^3 q^n}{1-q^n} = 1 + 240 \sum_{n=1}^\infty \sigma_3^{\phantom.}(n) q^n = 1 + 240 q + 2160 q^2 + 6720 q^3 + \cdots,
$$
$$
E_6 = 1 - 504 \sum_{n=1}^\infty \frac{n^5 q^n}{1-q^n} = 1 - 504 \sum_{n=1}^\infty \sigma_5^{\phantom.}(n) q^n = 1 - 504 q - 16632 q^2 - 122976 q^3 - \cdots,
$$
where $\sigma_k^{\phantom.}(n) = \sum_{d|n} d^k$, the sum of the $k$-th powers of the divisors of $n$.  Then $j = E_4^3 / \Delta$, where
$$
\Delta = \sum_{n=1}^\infty \tau(n) q^n = q - 24 q^2 + 252 q^3 \cdots
$$
may be computed either as $12^{-3} (E_4^3 - E_6^2)$ or from the product formula $\Delta = \eta^{24} = q \prod_{n=1}^\infty (1-q^n)^{24}$.
[Warning: Serre uses $E_k$, not the usual $E_{2k}$, for the normalized Eisenstein series of weight $2k$, so the above $E_4$ and $E_6$ are Serre's $E_2$ and $E_3$.  EDIT and I see that Álvaro Lozano-Robledo uses yet another convention, $g_2(\tau)$ and $g_3(\tau)$. TIDE]
This is also implemented in various computer algebra packages; e.g. in gp the command
ellj(q)

returns the power series $q^{-1} + 744 + 196884 q + \cdots$ out to whatever's the current series precision (by default $16$, so it ends $25497827389410525184 q^{14} + O(q^{15})$ because the next term would be of order $q^{16}$ times the leading term $q^{-1}$).
A: In addition to the power-series methods using Eisenstein series and $\Delta$, and the modular equation methods using, e.g., $h_5$ given in the other answers and comments, there are transcendental methods that allow you to compute individual coefficients numerically without working with the whole series.
H. Petersson gave an explicit formula for each $q$-expansion coefficient of $j$ in: Über die Entwicklungskoeffizienten der automorphen Formen, Acta Math. 58 (1932), 169–215.  Rademacher gave an independent derivation of the formula using his enhancement of the Hardy-Ramanujan circle method.  If we write $j(\tau) = \sum_{n \geq -1} c_n q^n$ with $c_{-1} = 1$, then for $n \geq 1$:
$$c_n = \frac{2\pi}{\sqrt{n}} \sum_{k = 1}^\infty \frac{A_k(n)}{k} I_1 \left( \frac{4\pi \sqrt{n}}{k} \right), \qquad \text{where} \quad A_k(n) = \sum_{ab \equiv -1 \pmod k} e^{\frac{2\pi i}{k}(na+b)},$$
and $I_1$ is a modified Bessel function of the first kind.  While this is an infinite sum, in practice you only need to add the first few terms because you already know the coefficients are integers.
More generally, if you want coefficients of any (weakly holomorphic) modular form of non-positive weight, you can apply the circle method, and you need only input the principal parts of the expansions at cusps.  The convergence seems to be faster for forms of "more negative" weight, e.g., the weight -12 form $\Delta^{-1}$ converges spectacularly quickly.  If you are just looking for expansions of Hauptmoduln such as the Monstrous Moonshine functions $T_g$, there is a formula for any genus zero group (that looks eerily similar to the formula for $j$) at the end of Venkov's 1983 ICM lecture.
