A radial basis interpolation function is described as:
$ f(\textbf{x})=\sum_{k=1}^N c_k \phi(\lVert \textbf{x}-\textbf{x}_k \rVert_2), \ \textbf{x}\in\mathbb{R}^s $
where $\textbf{x}_k$ are the $N$ scattered points and $c_k$ are the coefficients of the function obtained by solving the linear system:
$ \begin{bmatrix} f(\textbf{x}_1) \\ f(\textbf{x}_2) \\ \vdots \\ f(\textbf{x}_N) \end{bmatrix} = \begin{bmatrix} \varphi\left(\left\lvert \frac{1}{\lVert\textbf{x}_1\rVert} - \frac{1}{\lVert\textbf{x}_{1}\rVert} \right\rvert\right) & \varphi\left(\left\lvert \frac{1}{\lVert\textbf{x}_1\rVert} - \frac{1}{\lVert\textbf{x}_{2}\rVert} \right\rvert\right) & \ldots & \varphi\left(\left\lvert \frac{1}{\lVert\textbf{x}_1\rVert} - \frac{1}{\lVert\textbf{x}_{N}\rVert} \right\rvert\right) \\ \varphi\left(\left\lvert \frac{1}{\lVert\textbf{x}_2\rVert} - \frac{1}{\lVert\textbf{x}_{1}\rVert} \right\rvert\right) & \varphi\left(\left\lvert \frac{1}{\lVert\textbf{x}_2\rVert} - \frac{1}{\lVert\textbf{x}_{2}\rVert} \right\rvert\right) & \ldots & \varphi\left(\left\lvert \frac{1}{\lVert\textbf{x}_2\rVert} - \frac{1}{\lVert\textbf{x}_{N}\rVert} \right\rvert\right) \\ \vdots & \vdots & \ddots & \vdots \\ \varphi\left(\left\lvert \frac{1}{\lVert\textbf{x}_N\rVert} - \frac{1}{\lVert\textbf{x}_{1}\rVert} \right\rvert\right) & \varphi\left(\left\lvert \frac{1}{\lVert\textbf{x}_N\rVert} - \frac{1}{\lVert\textbf{x}_{2}\rVert} \right\rvert\right) & \ldots & \varphi\left(\left\lvert \frac{1}{\lVert\textbf{x}_N\rVert} - \frac{1}{\lVert\textbf{x}_{N}\rVert} \right\rvert\right) \\ \end{bmatrix} \begin{bmatrix} c_1 \\ c2 \\ \vdots \\ c_N \end{bmatrix} $
The radial basis function $\phi$ can be a Gaussian, an inverse multiquadric, etc. For a Gaussian, we get:
$ \phi(x) = e^{-(\epsilon x)^2} $
for some free parameter $\epsilon$
How do you find the derivative of the function and is it well-defined on all values or do we get into zero-denominator situations? for example, what is $\frac{df(x)}{dy}$ if $\textbf{x}\in\mathbb{R}^2$?
The coefficients $c_k$ are independent of $\mathbf x$. Therefore, the partial derivatives of $f$ are $$\frac{\partial}{\partial x_j}f(\textbf{x})=\sum_{k=1}^N c_k \frac{\partial}{\partial x_j}\phi(\lVert \textbf{x}-\textbf{x}_k \rVert_2)$$ By the chain rule, $$\frac{\partial}{\partial x_j}\phi(\lVert \textbf{x}-\textbf{x}_k \rVert_2) = \frac{\langle \mathbf e_j, \mathbf x-\mathbf x_k\rangle }{\|\mathbf x-\mathbf x_k\|_2}\phi'(\lVert \textbf{x}-\textbf{x}_k \rVert_2) $$ where $\mathbf e_j$ is the $j$th standard basis vector. That is, $\langle \mathbf e_j, \mathbf x-\mathbf x_k\rangle$ is just the $j$th component of $\mathbf x-\mathbf x_k$.
You can simplify this a bit by writing $\phi(x)=\psi(x^2)$. Then the argument of $\psi$ is squared norm, which is more pleasant to differentiate, the result being $$\frac{\partial}{\partial x_j}f(\textbf{x})=2 \sum_{k=1}^N c_k \langle \mathbf e_j, \mathbf x-\mathbf x_k\rangle \psi'(\lVert \textbf{x}-\textbf{x}_k \rVert_2)$$
In any event, $f$ has derivative everywhere, as long $f$ itself is defined.