The definition of an entangled state is one that cannot be decomposed into an arbitrary Kronecker product of individual state vectors corresponding to individual quantum systems. Mathematically speaking, an un-entangled state $|\psi\rangle$ is one that can be written in the following way:
$$|\psi\rangle \ = \ \displaystyle\bigotimes_n |\psi_n\rangle$$
Where $|\psi_n\rangle$ is the state of the $n$-th qudit, whatever you want. In the case of an $D$-level discrete quantum system, we can write each $|\psi_n\rangle$ as a normalized linear combination of $D$ orthogonal basis vectors:
$$|\psi_n\rangle \ = \ \displaystyle\sum_{d} c^n_d |d\rangle \ \ \ \text{with} \ \ \ \displaystyle\sum_d |c^n_d|^2 \ = \ 1$$
In the case of a qubit, $D \ = \ 2$ and we have $|0\rangle$ and $|1\rangle$ as our basis states. We thus expand $|\psi\rangle$:
$$|\psi\rangle \ = \ \displaystyle\bigotimes_n \displaystyle\sum_{d} c^n_d |d\rangle$$
When we expand our tensor product, we will get a linear combination of $D^N$ different state vectors, each with some unknown coefficient. In general, for the $j$-th term of the sum, our coefficient $c_j$ is given by a product of $D^N$ different values of $c^n_d$ (I'm not going to do the explicit expansion for the $j$-th term in the case of a $d$-level quantum system, I don't even want to do the binomial expansion for a $2$-level system):
$$c_j \ = \ c^{n_{j_1}}_{d_{j_1}} \ c^{n_{j_2}}_{d_{j_2}} \ ... \ c^{n_{j_{D^N}}}_{d_{j_{D^N}}}$$
But we already know $c_j$, as we know $|\psi\rangle$, and $c_j$ is given by:
$$c_j \ = \ \langle j | \psi \rangle$$
Where $|j\rangle$ is the $j$-th basis vector in the expansion. Thus, all we have to do is solve the system of $D^N$ equations for $D^N$ unknowns. If we can in fact do this, then we know that we have an un-coupled (un-entangled) state. If we can't we know our state is entangled!
If this was slightly confusing an abstract, let's work through a basic example, and prove that the maximally entangled qubit Bell state:
$$|\psi\rangle \ = \ \frac{|00\rangle \ + \ |11\rangle}{\sqrt{2}}$$
is in fact entangled using the outlined method. Our goal is to write $|\psi\rangle$ as a tensor product of two individual state vectors:
$$|\psi\rangle \ = \ \frac{|00\rangle \ + \ |11\rangle}{\sqrt{2}} \ = \ (a|0\rangle \ + \ b|1\rangle) \ \otimes \ (c|0\rangle \ + \ d|1\rangle)$$
$$\Rightarrow \ (a|0\rangle \ + \ b|1\rangle) \ \otimes \ (c|0\rangle \ + \ d|1\rangle) \ = \ ac |00\rangle \ + \ ad |01\rangle \ + \ bc|10\rangle \ + \ bd|11\rangle$$
Now, since each of our basis vectors are orthogonal, we can equate the components to get a system of equations:
$$ac \ = \ \frac{1}{\sqrt{2}}$$
$$ad \ = \ 0$$
$$bc \ = \ 0$$
$$bd \ = \ \frac{1}{\sqrt{2}}$$
Rearranging, we get:
$$d \ = \ \frac{1}{\sqrt{2} b} \ \ \ \text{and} \ \ \ a \ = \ \frac{1}{\sqrt{2} c}$$
So we get:
$$ad \ = \ \frac{1}{2bc} \ = \ 0$$
Which is a contradiction (this is not possible). Thus, we cannot express the Bell state as a tensor product of individual state vectors. It follows that the Bell state is in fact entangled.
If you are into some really hardcore stuff, we can do something similar with a continuous-variable quantum system (qumodes, for instance). In this case, we no longer expand out our state $|\psi_n\rangle$ as a discrete sum of basis states, but as a continuous integral over a continuum:
$$|\psi_n\rangle \ = \ \displaystyle\int \psi_n(x) |x\rangle dx$$
We then have:
$$|\psi\rangle \ = \ \displaystyle\bigotimes_n \displaystyle\int \psi_n(x) |x\rangle dx \ = \ \displaystyle\int \ ... \ \displaystyle\int \displaystyle\bigotimes_n \psi_n(x_n) |x_n\rangle \ dx_1 \ ... \ dx_n$$
In general, $|\psi\rangle$ is given by:
$$|\psi\rangle \ = \ \displaystyle\int \ ... \ \displaystyle\int \psi(x_1, \ ..., \ x_n) \displaystyle\bigotimes_n |x_n\rangle \ dx_1 \ ... \ dx_n$$
And, analogous to the discrete case, we know the function $\psi(x_1, \ ..., \ x_n)$. Thus, in the un-entangled case, we have:
$$\displaystyle\prod_n \psi_n(x_n) \ = \ \psi(x_1, \ ..., \ x_n)$$
So instead of finding discrete components, our goal is to find a collection of functions $\psi_j(x_j)$, such that the product of all these functions gives us the general $\psi$. If we can't do this, then our state is entangled. As an example, consider the function:
$$\psi(x, \ y) \ = \ x \ + \ y \ + \ xy$$
We want to find $\psi_1$ and $\psi_2$ such that:
$$\psi_1(x) \psi_2(y) \ = \ x \ + \ y \ + \ xy$$
We have:
$$\psi_2(y) \ = \ \frac{x \ + \ y \ + \ xy}{\psi_1(x)} \ = \ \frac{x}{\psi_1(x)} \ + \ \frac{y}{\psi_1(x)} \ + \ \frac{xy}{\psi_1(x)}$$
Now, $\psi_2(y)$ is a function of only $y$, we can't have $x$ present in any of the terms. We quickly realize that this is impossible, as the only way to get rid of the $x$ in the first term is if $\psi_1(x) \ = \ ax$, but this would mean that the second term contains $x$. This same logic can be applied to all of the terms. In addition, cancellation is obviously of no help here. Thus, the state given as:
$$|\psi\rangle \ = \ \displaystyle\int \displaystyle\int (x \ + \ y \ + \ xy) |x, \ y\rangle dx dy$$
must be entangled.
Edit
After writing this answer, I came across this question, which is very similar, and already has answers. Despite this, I don't think this question should be closed as my answer discusses a more general approach to this problem.
