How to compute a Differential Equation With Dot Product. I have the following differential equation, that I would like to plug into a numerical integrator.
$$
\frac{d \mathbf r}{dt} \cdot \hat{n}(\mathbf r) = 0
$$
Where $\mathbf r(t)$, is the trajectory of a particle and $\hat n$ is a normal vector to a surface. So this equation says that the velocity of the particle must be perpendicular to the normal of the surface. It would be useful to have $d\mathbf r / dt$ alone. To achieve this we impose other restrictions like:
$$
\left|\left| \frac{d \mathbf r}{dt} \right|\right| = 1
$$
We get two equations:
$$
\begin{pmatrix}
n_x & n_y & n_z\\
u_x & u_y & u_z\\
1 & 1 & 1
\end{pmatrix}
\begin{pmatrix}
u_x \\
u_y \\
u_z \\
\end{pmatrix} =
\begin{pmatrix}
0 \\
1 \\
u_x + u_y + u_z \\
\end{pmatrix}
$$
And we let the last equation be whatever we want, obtaining the following result:
$$
\mathbf u(t) \equiv \frac{d \mathbf r}{dt}
$$
$$
\mathbf u(t) =\begin{pmatrix}
u_z - u_y & n_y - n_z & u_y n_z - u_z n_y\\
u_x - u_z & n_z - n_x & u_z n_x - u_x n_z\\
u_y - u_x & n_x- n_y & u_x n_y - u_y n_x
\end{pmatrix} \begin{pmatrix}
0 \\
1 \\
u_x + u_y + u_z \\
\end{pmatrix} = 
\begin{pmatrix}
n_y -n_z + (u_x + u_y +u_z) (u_yn_z-u_zn_y) \\
n_z -n_x + (u_x + u_y +u_z) (u_yn_z-u_zn_y) \\
n_x -n_y + (u_x + u_y +u_z) (u_yn_z-u_zn_y) \\
\end{pmatrix}
$$
$$
\mathbf u(t) = (\mathbf i + \mathbf j + \mathbf k) \times \hat n + ((\mathbf i + \mathbf j + \mathbf k) \cdot \mathbf u(t)) (\mathbf u(t) \times\hat n) 
\\
(\mathbf i + \mathbf j + \mathbf k) \equiv \mathbf e \\
\mathbf e \cdot \mathbf u(t) \equiv k(t) \\
\mathbf u(t) = (\mathbf e  + k(t)\mathbf u(t)) \times \hat n
\\
$$
This last equation is equivalent (I suppose), to the first one I proposed.
But how would I plug this into a program. What I'm seeking for is an equation of the form,
$$
\mathbf u_{i + 1} = F(\mathbf u_i, \mathbf r_i)
$$
that relates the previous $\mathbf u(t)$ with the new one, but I do not have $\mathbf u(t)$ alone (basically because I can't). How I'm supposed to do this?

There is not enough data to compute the trajectory. The equation you present just says that the velocity will always be tangent to the surface where the motion occurs. This must be supplemented with the actual equations of motion. A natural choice would be
$$
m\, r''(t) = (0, 0, -g).
$$
