I don't think it's possible to tell if a function is continuous/differentiable/smooth without knowing the function definition, i.e. its code. I don't even know how these concepts are defined in finite precision math. After all, the set of floating point numbers in a computer is not continuous.
On the other hand, it is easy to calculate the Jacobian of a function numerically. A Jacobian is a bunch of partial derivatives, and you can calculate them with finite differences. For example:
$$\mathbf{F}(\mathbf{x}) = u(x, y, z) \mathbf{i} + v(x, y, z)\mathbf{j} + w(x, y, z)\mathbf{k}$$
$$J = \begin{bmatrix}
\dfrac{\partial u}{\partial x} & \dfrac{\partial u}{\partial y} & \dfrac{\partial u}{\partial z}
\\
\dfrac{\partial v}{\partial x} & \dfrac{\partial v}{\partial y} & \dfrac{\partial v}{\partial z}
\\
\dfrac{\partial w}{\partial x} & \dfrac{\partial w}{\partial y} & \dfrac{\partial w}{\partial z}
\end{bmatrix}$$
$$\frac{\partial u}{\partial x} \approx \frac{u(x+h, y, z) - u(x-h, y, z)}{2h}$$
$h$ is a small number like $10^{-3}$. If it is too big, you will get big truncation errors. If it is too small, you will get big roundoff errors. Usually it must be calibrated by trial and error until you get good results.
Often what you need is the Jacobian determinant, which can be calculated from the above matrix.
When using finite differences, you don't need to know the function definition. You just need to call it a couple of times with different parameters.