A more elegant way of computing this Wronskian? As I was working on my differential equation homework this week I came across this problem:
Let $y^{(4)} + 16y=0$. Compute the Wronskian of four linearly independent solutions.
It's rather straightforward to find four such solutions solutions:
$\phi_1(x)= e^{\sqrt{2}x} \cos{\sqrt{2}x}$,
$\phi_2(x)= e^{\sqrt{2}x} \sin{\sqrt{2}x}$,
$\phi_3(x)= e^{-\sqrt{2}x} \cos{\sqrt{2}x}$,
$\phi_4(x)= e^{-\sqrt{2}x} \sin{\sqrt{2}x}.$
And from there computing the Wronskian (it is $256$) can be accomplished by trudging through the computation. 
However, since the solutions have such a nice symmetry on the complex unit circle, is there an easier way of coming to the solution? I tried to do this problem at first by avoiding taking the determinant, but ended spending more time than I would have anyway in doing so. I want to believe there is an easier way to do this problem, and I feel like there is a trick that I'm not seeing.
 A: 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.
$$
A: 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).$$
