How to compute the Helmholtz decomposition of 2D and 3D vector fields?

I have a sample of 100 particles with their 3D positions and velocities. I want to decompose the velocity vector field to its curl free and divergence free components.

Is it a feasible to do so for the sampled data using Python, Matlab or Mathematica?

Certainly, you can try to find the Helmholtz decomposition on your sampled data, and find your irrotational and solenoidal components. However, there are certain requirements on your original vector field you started with. In general, it is that the vector field you are trying to decompose has to be sufficiently smooth and decay rapidly. Commonly, that corresponds to the vector field $\mathbf F$ being twice continuously differentiable on some bounded domain $V$, and $\mathbf F$ decay faster than $1/r$, leading to $$\mathbf F=-\nabla\Phi+\nabla\times\mathbf A$$ decomposition with scalar potential $\Phi$ (for the irrotational component) and vector potential $\mathbf A$ (for the solenoidal component).