As with all derivations, it depends on what you want to treat as fundamental. Typically we would derive Coulomb's law from the Maxwell equations, so we're trying to solve
$$\nabla\cdot \mathbf{E} = -\nabla^2 \varphi = q\delta(\mathbf{x})/\epsilon_0\qquad (1)$$
In $n$ spatial dimensions and in Cartesian coordinates $(x_1,\ldots,x^n)$, this becomes
$$\sum_{k=1}^n \frac{\partial^2}{\partial x_k^2} \varphi = -\frac{q}{\epsilon_0}\delta(\mathbf x)\qquad(2)$$
Because this problem has spherical symmetry, we can move to hyperspherical coordinates. If we do so, we will find that$^\dagger$
$$\frac{1}{r^{n-1}}\frac{\partial }{\partial r}\left(r^{n-1} \frac{\partial\varphi}{\partial r}\right) = -\frac{q}{\epsilon_0} \delta(r)\qquad (3)$$
Away from $r=0$, we would therefore have that $$\frac{\partial}{\partial r}\left(r^{n-1} \frac{\partial \varphi}{\partial r}\right)=0 \implies r^{n-1} \frac{\partial \varphi}{\partial r} = c$$ for some constant $c$, and therefore that $\varphi = c\ r^{2-n}+d$ (unless $n=2$, in which case we'd have a logarithm). The constant $d$ can be set to zero by demanding that the potential vanish at infinity (this is an arbitrary choice, but a convenient one). The constant $c$ can be determined by using the divergence theorem to integrate $(1)$ over a hypersphere of radius $R$. Because of the spherical symmetry, the left hand side would be the surface area of the $(n-1)$-sphere of radius $R$ times $\varphi'(R)$:
$$-\frac{2\pi^{n/2}}{\Gamma(n/2)}R^{n-1} \varphi'(R)=\left(\frac{2\pi^{n/2}(n-2)}{\Gamma(n/2)}\right) c$$
while the right hand side is simply equal to $q/\epsilon_0$ because of the delta function. As a result,
$$\varphi(r) = \frac{\Gamma(n/2)}{2(n-2)\pi^{n/2}\epsilon_0} \frac{q}{r^{n-2}}\qquad (4)$$
In $n=3$ dimensions, we have $\Gamma(3/2)=\sqrt{\pi}/2$ so this reduces to the familiar case
$$\varphi^{(3)}(r) = \frac{1}{4\pi\epsilon_0} \frac{q}{r} \implies \mathbf{E}^{(3)}(r) = \frac{1}{4\pi\epsilon_0} \frac{q}{r^2}\hat r$$
In 4-dimensions, $\Gamma(2)=1$ so we would have
$$\varphi^{(4)}(r) = \frac{1}{4\pi^2 \epsilon_0} \frac{q}{r^2} \implies \mathbf{E}^{(4)}(r) = \frac{1}{2\pi^2 \epsilon_0} \frac{q}{r^3} \hat r$$
In the other direction, for $n=1$ we have $\Gamma(1/2)=\sqrt{\pi}$ and so
$$\varphi^{(1)}(r) = -\frac{1}{2\epsilon_0} q r \implies \underbrace{\mathbf{E}^{(1)}(r)=\frac{1}{2\epsilon_0} q \hat r}_{\text{constant}}$$
So intuitively the equation should not change?
The problem is that $\nabla^2$ changes in higher dimensions, so if you reuse the familiar form of Coulomb's law then it will not obey the Maxwell equations. Assuming you'd like to treat the latter as more fundamental, we need to use Gauss' law to find the more general form of Coulomb's law.
Is Gauss law applicable only in 3 dimensions or valid for any dimension, because the Gauss divergence theorem is only for 3 dimensions?
The divergence theorem holds in an arbitrary number of dimensions. If we assume that Gauss' law holds in an arbitrary number of dimensions, then we find Coulomb's law as I did above. Of course, Gauss' law is a physical statement, not a purely mathematical one, so there's no way to mathematically prove that it holds for all dimensions.
$^\dagger$This expression should not be taken too literally, as the delta function at the origin has some pathological issues in spherical coordinates. The spirit of this equation is that we will find the solution for $r\neq 0$, and obtain the remaining undetermined constant by integrating $(1)$.