How to create a vector field whose Curl and Divergence are zero at any point?

Take any harmonic function $\phi$, and set your vector field to be $v = \nabla \phi$, the gradient.

Edit On any simply connected domain (if your space is just the entire Euclidean space, for example), any curl free vector field can be written as the gradient of a function. You can see this explicitly by verifying that given a curl free vector field $v$, the function $\phi(x) = \int_{C} v(s)\cdot ds$ where $C$ is a curve starting at the origin and ending at the point $x$ is independent of the choice of the curve $C$, and hence $\phi$ is well-defined. (This is generally treated in any introductory textbook on vector calculus, and is a special instance of Stokes' theorem.) Once you have that, plugging $v = \nabla \phi$ into the divergence equation you automatically get that $\nabla^2\phi = 0$.

There is a theorem (which apriori has nothing to do with Stokes' theorem) which is called Poincarais' Lemma. The version that applies here states that in $U$ - a sufficiently convex subset of (or whole) $R^3$.

a) for every divergence-free vector field V there exists another field A such that $\nabla \times A = V$.

b) for every curl-free vector field V there exists scalar field $\phi$ such that $\nabla \phi = V$.

Consult textbooks if interested in definition of 'sufficiently convex'.

One can use one of those statements to simplify our search - because using this theorem reduces our requirements from two ($\nabla \times V = 0, \nabla \cdot V = 0$) to one. For example, let us use b). Then, our search reduces to finding a solution to Laplace's equation $\nabla \cdot \nabla \phi = 0$ and the second requirement ($\nabla \times V = 0$) will be automatically satisfied. Furthermore, Poincarais' Lemma makes us certain that this is a way of generating ALL such fields.

Now, how to solve Laplace equation?

It appears that well-posedness (existence and uniqueness of a solution) of this equation requires additional information - the boundary data. More specifically, we have to supply smooth values of field $\phi$ on the boundary $\partial U$ of the set for which Poincarais' Lemma works.

To understand this, let us see some examples:

1) $U = R^3$ (whole space). Laplace's equation reads $$\left(\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}+\frac{\partial^2}{\partial z^2}\right) \phi(x,y,z) = 0$$ boundary values - $\phi = 0$ at all infinities. Then we immediately see, that we can find a solution which satisfies those boundary conditions - namely $\phi = 0$ everywhere! Because solution has to be unique, we conclude that the field found in this case is $\nabla \phi = \vec{0} $ which is of course true - the 0 vector has no divergence and curl!

2) $U = R^3$ but this time we demand $\phi \to \infty$ as $x \to \infty$ and $\phi \to -\infty$ as $x \to -\infty$. We do not say a thing about other boundary layers so the boundary value problem is singular in this case - fortunately what we lose is just uniqueness. For example, we can take $\phi(x,y,z)=x$ as well as $\phi(x,y,z)=3x + y$.

The corresponding fields will be $(1,0,0)$ and $(3,1,0)$.

To classify solutions we have to analyze all the possible boundary values. Yes, it is a pain. But for most real-world fields which decay at infinity the answer is simply - there are no such fields unless they are $\vec{0}$ - case 1)

When $U$ is finite, one has to impose boundary conditions (which is easier than on unbounded domain) which most often come directly from the physical problem in hand. Then, one has to adapt coordinates suitable for such boundary structure (spherical for sphere, cartesian for box, etc.) and solve BVP for Laplace equation - which is different every time - no general solution. However, most textbooks on electrodynamics are full of examples.

The most interesting part (and least seen in physics) is when $U$ is so chosen that the assumptions of Poincarais Lemma are not satisfied. But for this, consult mathematical exchange, please :)