For the circular paraboloid $$\Sigma : z=x^2+y^2,$$ I was able to find a Gauss map $N : \Sigma\to S^2$ explicitly as $$N((x, y, z))=\left(-\dfrac{2x}{\sqrt{4z+1}}, -\dfrac{2y}{\sqrt{4z+1}}, \dfrac{1}{\sqrt{4z+1}}\right).$$
How can I compute the shape operator of $\Sigma$ at some point, say $p=(1, 1, 2)$?
Hmm, here is how I'd do it. So first let us find the tangent vectors from the positiosn vectors $R(x,y,z) = (x,y,x^2+y^2)$:
$$ \tau_1= \partial_x R = (1,0,2x)$$
$$ \tau_2= \partial_y R = (0,1,2y)$$
Now to know a matrix, it's sufficient to know it's action on the basis vectors. We know that $-\nabla_{\tau_1} n = S(\tau_1)$ and $-\nabla_{\tau_2} n = S_{(\tau_2)}$. Also note that:
$$ \nabla_{v} n = \nabla_{v} (n_x e_x + n_y e_y +n_ze_z) = (\nabla_v n_x)e_x + (\nabla_v n_y) e_y + (\nabla_v n_z) e_z$$
The $\nabla_v$ doesn't effect the basis because the cartesian basis is the same everywhere in space.
Putting everything together, you should have your answer.