how do I compute the fundmental domain for a congruence subgroup of $SL(2, \mathbb{Z})$
This region is important because the theta function $\theta(z) = q^{n^2}$ is invariant under two transformations $z \mapsto z + 1$ and $z \to - \frac{1}{4z}$, and these two generate the congruence subgroup $\Gamma_0(4)$.
I know there is an algorithm for any Fuchsian group. Here is an erroneous comutation of the fundamental domain. It has too many cusps.
I found a blog with a picture of what could be the fundamental domain of $\Gamma_0(2)$. Therefore the domain for $\Gamma_0(4)$ could be related. I know that $[SL(2,\mathbb{Z}): \Gamma_0(2)]= 3$ and $[SL(2,\mathbb{Z}): \Gamma_0(4)]= 6$ but I can't figure out which two copies to join. And ther's no derivation for this particular case.
As suggested in the comments, we can obtain a fundamental domain for $\Gamma_0(4)$ by taking a fundamental domain for $\newcommand{\SL}{\operatorname{SL}} \newcommand{\Z}{\mathbb{Z}} \SL_2(\Z)$ and translating by coset representatives for $\SL_2(\Z)/\Gamma_0(4)$.
First let's recover the fundamental domain for $\Gamma_0(2)$ that you linked. We begin with the fundamental domain $D_1$ for $\SL_2(\Z)$ below.
$\hspace 6cm$
Taking the coset representatives \begin{align*} \left(\begin{array}{rr} 1 & 0 \\ 0 & 1 \end{array}\right), \left(\begin{array}{rr} 0 & 1 \\ -1 & 1 \end{array}\right), \left(\begin{array}{rr} -1 & 1 \\ -1 & 0 \end{array}\right) \end{align*} for $\SL_2(\Z)/\Gamma_0(2)$ and translating $D_1$ by each element, we obtain the fundamental domain $D_2$ for $\Gamma_0(2)$ below.
$\hspace 6cm$
To obtain a fundamental domain for $\Gamma_0(4)$, we choose a representative $A$ for the nontrivial coset of $\Gamma_0(2)/\Gamma_0(4)$ and then form $D_2 \cup (A \cdot D_2)$.
To get a fundamental domain similar to the one in the picture you drew, we want an element that acts on $\frac{1}{2} + \frac{1}{2} i$ in the same way that $S = \begin{pmatrix} 0 & -1\\ 1 & 0\end{pmatrix}$ acts on $i$. So we conjugate $S$ by the element $\left(\begin{array}{rr} 0 & 1\\ -1 & 1\end{array}\right)$ that takes $i$ to $\frac{1}{2} + \frac{1}{2} i$: \begin{align*} A := \left(\begin{array}{rr} 0 & 1\\ -1 & 1\end{array}\right) \begin{pmatrix} 0 & -1\\ 1 & 0\end{pmatrix} \left(\begin{array}{rr} 0 & 1\\ -1 & 1\end{array}\right)^{-1} = \left(\begin{array}{rr} 1 & -1\\ 2 & -1\end{array}\right) \, . \end{align*} Note that $A \in \Gamma_0(2) \setminus \Gamma_0(4)$, as required. Multiplying each of the coset representatives for $\Gamma_0(2)$ by $A$, we obtain the coset representatives \begin{align*} \left(\begin{array}{rr} 1 & 0 \\ 0 & 1 \end{array}\right), \left(\begin{array}{rr} 0 & 1 \\ -1 & 1 \end{array}\right), \left(\begin{array}{rr} -1 & 1 \\ -1 & 0 \end{array}\right), \left(\begin{array}{rr} 1 & -1 \\ 2 & -1 \end{array}\right), \left(\begin{array}{rr} 1 & 0 \\ 1 & 1 \end{array}\right), \left(\begin{array}{rr} 0 & 1 \\ -1 & 2 \end{array}\right) \end{align*} with the corresponding fundamental domain below.
$\hspace 6cm$
As a note, if you only need the area of a fundamental domain for $\Gamma_0(4)$, we don't actually have to construct the fundamental domain. Since $[\SL_2(\Z) : \Gamma_0(4)] = 6$, then $$ \newcommand{\vol}{\operatorname{vol}} \vol(\Gamma_0(4)) = [\SL_2(\Z) : \Gamma_0(4)] \vol(\SL_2(\Z)) = 6 \cdot \frac{\pi}{3} = 2\pi \, . $$
(All these images were made with SageMath. Here's a link to a SageMathCell with the commands.)