How to create a vector field whose Curl and Divergence are zero at any point? What is the mathematical procedure to derive a vector field whose curl and divergence are zero at any point at any time?
Edit: Please explain it by solving the differential equations of curl and divergence.
 A: 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 :)
A: 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$. 
A: A nice example in physics is the standard expression for a dipole field, either magnetic or electric. With the exception of the origin, where the dipole is located, div and curl are zero, and so must also satisfy the Laplace equation. But, wait! del$^2$$\left(\frac{1}{r}\right)$ is $-4\pi$ (Dirac $\delta(r)$) ...ah grass hopper...
A: A differentiable vector field $F=(u,v)$ defined on an open set $\Omega\subset\mathbb{R^2}$ is curl and divergence free if and only if $f=u-iv$ is a holomorphic function (easy consequence of Cauchy-Riemann equations).  
Then, if you want to obtain a differentiable curl and divergence free vector field on an open set $\Omega\subset\mathbb{R^2}$, simply take the conjugate of any holomorphic function defined on $\Omega$.
A: Cleaning up for us lazy physicists.
1) First solution:
We want: $$\nabla\cdot {\vec{B}}=0$$$$\nabla\times {\vec{B}}=0$$
From the above mentioned theorem of Poincare:
$$\nabla\cdot {\vec{B}}=0 \implies \vec{B}=\nabla \times\vec{A}$$
By the identity:
$$\nabla \times(\nabla \times\vec{A})=\nabla(\nabla \cdot \vec{A})-\nabla^2\vec{A}$$
So the divergence and curl are both zero if we have $\vec{B}=\nabla \times\vec{A}$ and $\nabla(\nabla \cdot \vec{A})-\nabla^2\vec{A}=0$
2) Seconds solution.
Again by Poincare as mentioned above:
$\nabla \times \vec{B}=0 \implies \vec{B}=\nabla \phi$
Then we have $\nabla \cdot \vec{B}=0$ if:  $$\nabla^2 \phi=0$$
This is Laplace's equation. Probably the most important you'd encounter in undergrad physics. 
Legendre Polynomials are among the most important solutions to this equation. They represent solutions to Maxwell's Equations for static charges and currents as well as angular momentum in quantum mechanics. 
