The space of cusp forms on $\mathbb{H}/\Gamma_0(4)$ is finite dimensional and spanned by the Poincare series for $m \in \mathbb{Z}$:
$$ P_m \big(z; \Gamma_0(4) \big) = \sum_{\tau \in \Gamma_\infty \backslash \Gamma_0(4) } j(\tau, z)^{-2k} \, e\,\big(\,m \tau z\,\big) = \sum_{\tau \in \Gamma_\infty \backslash \Gamma_0(4) } (\pm i ) \, \frac{1}{ | \, cz + d \, |^k } \, e\,\big(\,m \tau z\,\big)$$
The theta multiplier $j(\tau, z)$ is describing how theta functions transform: $$ j(\tau, z) = \theta(\tau z)/\theta(z) = \varepsilon_d^{-1} \left( \frac{c}{d}\right) \big( cz + d \big)^{1/2} \quad\text{ with }\quad\varepsilon = \left\{ \begin{array}{cc} 1 & \text{if }d=4k+1\\ i & \text{if }d=4k-1\end{array} \right. $$
One question I have is how can the space $S_k\big(\Gamma_0(4)\big)$ be finite dimensional, when there are infinitely many Poincare series, $m \in \mathbb{Z}$. There must be infinitely many relations between the different Poincare series $P_m$.
$$ \Gamma_\infty = \left\{ \left( \begin{array}{cc} 1 & x \\ 0 & 1 \end{array} \right) : x \in \mathbb{Z} \right\} $$
Before I continue, I need to know how to write down representatives of $\Gamma_\infty \backslash \Gamma_0(4)$. Are these the same as the "cusps" of $\Gamma_0(4)$? I was able to find related questions for that:
To add to the confusion, I'd read that spaces of modular forms split into the Eistenstein and cusp parts:
$$ M_k\big(\Gamma_0(4)\big)=E_k\big(\Gamma_0(4)\big)\oplus\text{S}_k\big(\Gamma_0(4)\big) $$
and here I can get the Eisenstein series, by setting $m = 0$. If instead I was looking for Eistensin series over $SL(2, \mathbb{Z})$ the answer is in many textbooks:
$$ E(z,s) = \sum_{\Gamma_\infty\backslash \Gamma } \mathrm{Im}(\gamma z)^s = \sum_{(c,d) = 1} \frac{y^s}{ \big|cz + d \big|^{2s}} $$
and we can do this because of a Lemma; we have the bijection:
\begin{eqnarray} \Gamma_\infty \backslash \Gamma_0 &\to& \big\{ (x,y)\in \mathbb{Z}^2 /\pm 1 : \mathrm{gcd}(x,y) = 1 \big\} \\ \Gamma_\infty\left( \begin{array}{cc} a & b \\ c & d\end{array} \right) &\mapsto & \pm (c,d) \end{eqnarray}
What is the corresponding Lemma for $\Gamma_\infty \backslash \Gamma_0(4)$ ? and what are the appropriate numbers $(c,d)$ to finish writing the Poincaré series on the first line?
Okay, this makes no sense. First of all, why are you using theta multipliers? This comes up in half-integral weight modular forms, not modular forms of integral weight.
Second of all, read basically anything by Iwaniec to understand what these Poincare series look like. The correct $j$ factor is $j(\gamma,z) = cz + d$ for $\gamma = \begin{pmatrix} a & b \\ c & d \end{pmatrix}$. Then the Poincare series is \[P_m(z) = \sum_{\gamma \in \Gamma_{\infty} \backslash \Gamma_0(q)} j(\gamma,z)^{-k} e(m\gamma z).\] The Bruhat decomposition implies that $\Gamma_{\infty} \backslash \Gamma_0(q)$ has a set of representatives consisting of the identity and $\begin{pmatrix} * & * \\ c & d \end{pmatrix}$ with $c \geq 1$, $c \equiv 0 \pmod{q}$, and $d \in (\mathbb{Z}/c\mathbb{Z})^{\times}$. This allows you to find that the Fourier coefficients of $P_m(z)$; see Lemma 14.2 of Iwaniec and Kowalski.
Yes, even though there are infinitely many Poincare series, they only span a finite-dimensional space, so there are tons of linear relations between them. In fact, you can calculate their inner product directly! (This is the right-hand side of the Petersson trace formula.)