".... the question is that the divergence equation is a first order PDE so by giving enough boundary conditions we should be able to determine the field right?"
Not so. As mentioned in the comments the answer to this question is essentially about the Helmholtz decomposition, but let's actually partly delve into a certain proof of this decomposition that shows very clearly, geometrically and intuitively what the problem is, at least for a large class of vector fields, namely, those that have Fourier transform, as discussed in my answer here.
Imagine a vector field's Fourier decomposition $\vec{F}(\vec{k})$, a function of the plane wave wavevector $\vec{k}$, i.e. we decompose a vector valued function $\vec{f}(\vec{r})$ of position $\vec{r}$ into a superposition of plane wave vector fields of the form $\vec{F}(\vec{k})\,\exp(i\,\vec{k}\cdot\vec{r})$.
Now, how do divergence and curl look in Fourier space? $\nabla\cdot \vec{f}$ has the Fourier transform $\vec{k}\cdot\vec{F}$ and $\nabla\times \vec{f}$ has the transform $\vec{k}\times\vec{F}$; you should be able to prove this fairly straightforwardly.
So now, ask your question in Fourier space terms. It is, "why can we determine the vector $\vec{F}$ from $\vec{k}\cdot\vec{F}$ alone?". It should be very clear that this can't be done; we need to know the components of $\vec{F}$ that are orthogonal to $\vec{k}$ and these can be assigned essentially independently, since the divergence of a vector field everywhere orthogonal to $\vec{k}$ vanishes.
In general, one can assign a smooth scalar field in Fourier space $g(\vec{k})$ and a second smooth vector field $\vec{H}(\vec{k})$ that is everywhere orthogonal to $\vec{k}$, but otherwise arbitrary. As I discuss in this answer here and here, the information $g(\vec{k})$ and $\vec{H}(\vec{k})$ are exactly the information to determine a vector field $\vec{F}$ such that:
$$g(\vec{k}) = \vec{k}\cdot\vec{F}(\vec{k})$$
$$\vec{H}(\vec{k})= \vec{k}\times \vec{F}(\vec{k})$$
So the answer to your question is essentially that the divergence condition tells you only the component of the vector field that is along the wavevector; the solenoidal part, that orthogonal to the wavevector, is missing (it has zero divergence) and can be assigned independently.