A more elegant way of computing this Wronskian?

We have the Liuville formula: $$ W(x)=W(x_0)\exp\left\{-\int_{x_0}^x\frac{a_1(x)}{a_0(x)}dx\right\} , $$ where $W(x)$ is the Wronskian, and $a_0(x)$ and $a_1(x)$ are $$ a_0(x)y^{(n)}+a_1(x)y^{n-1}+\ldots+a_n(x)y=0. $$ Your $a_1(x)$ is zero, hence $$ W(x)=W(x_0)=\rm const. $$ Since we are free to choose any constant, the answer is, e.g., $$ W(x)=1. $$

Added: One can determine the Wronskian up to a multiplicative constant using the equation coefficients, as explained in Artem's answer (my apologies: initially I misunderstood his point). Here this gives simply a constant function. If one is interested in determining this unknown constant explicitly for a particular basis of solutions, then the usual way (e.g. for special function ODEs) is to keep a sufficient number of terms in the asymptotics of solutions. For example here one could replace them by Taylor expansions truncated at 3rd order and compute the Wronskian dropping all non-constant terms. What follows is an alternative to that and to the direct approach.

---

One option is to compute instead the Wronskian of $\phi_{\pm}=\phi_{1} \pm i \phi_2$, $\psi_{\pm}=\phi_3\pm i\phi_4$. The main point is that $$\phi_{\pm}(x)=e^{(1\pm i)\sqrt{2} x},\qquad \psi_{\pm}(x)=e^{(-1\pm i)\sqrt{2} x},$$ and the derivatives of exponentials are easy to compute. One then finds \begin{align} &W(\phi_+,\phi_-,\psi_+,\psi_-)=-4W(\phi_1,\phi_2,\phi_3,\phi_4)=\\ &=\sqrt{2}\cdot\left(\sqrt{2}\right)^2\cdot\left(\sqrt{2}\right)^3 \operatorname{det}\left(\begin{array}{cccc} 1 & 1 & 1 & 1 \\ 1+i & 1-i & -1+i & -1-i \\ (1+i)^2 & (1-i)^2 & (-1+i)^2 & (-1-i)^2\\ (1+i)^3 & (1-i)^3 & (-1+i)^3 & (-1-i)^3 \end{array}\right)=\\ &=8\prod_{1\leq j<k\leq 4}\left(\lambda_k-\lambda_j\right)=\\ &=8\cdot (-2i)\cdot(-2)\cdot(-2-2i)\cdot(-2+2i)\cdot(-2)\cdot(-2i)=-1024. \end{align} In the 3rd line, we used the expression for the Vandermonde determinant with $$\left(\lambda_1,\lambda_2,\lambda_3,\lambda_4\right)=\left(1+i,1-i,-1+i,-1-i\right).$$