Thanks for the help, indeed it is similar to the Lyapunov equation.
Although for my case, the dimensions of $X$ and $B$ can be different from those of $G$ (that is: not square matrices, but of compatible dimensions).
But I can just solve the Lyapunov equation for each column of $B$ at a time and build my solution in that way (inside a loop). As such I don't even need to build the full matrix $G$ and can work with only the smaller ones $A$ and $B$.
\begin{eqnarray}
(I⊗A+A⊗I).X = B
\end{eqnarray}
$X,B \in \Re^{m \times n}$ and $(I⊗A+A⊗I) \in \Re^{m.m \times m.m}$
let $y = X(:,j)$ and $b = B(:,j)$ with $j = 1, 2 ... n$ (i.e. the columns of $X$ and $B$)
$v = (I⊗A+A⊗I).y \Leftrightarrow v = (I⊗A+A⊗I).vec(y)$
Knowing that: $(A⊗B).vec(X) = vec(BXA^T)$,
the previous expression can be transformed into $v = vec(AY + YA^T) $ or $V = AY + YA^T$
Thus we obtain the Lyapunov equation $AY + YA^T - V = 0$ or $AY + YA^T + \widetilde V = 0$
Now we set $\widetilde V = reshape(-b, m, m)$ (i.e. we change the vector $b$ into a matrix of dimensions $m \times m$) and solve the previous Lyapunov equation.
Finally $X(:,j) = vec(Y)$
Maybe there is a more straightforward way but this seems to work well for me now.