For Riemannian submanifolds, Eq. (5.41) (which gives the Riemannian Hessian in terms of the Euclidean Hessian + a correction term involving the Weingarten map) is informative theoretically but (as you saw in the Manopt code) we seldom use it "as is" to convert Euclidean derivatives into Riemannian derivatives. Rather, we work out the Riemannian Hessian "directly" and obtain a simplified formula "automatically". For example, see (7.29) for Stiefel. Of course, you would obtain the same formula via (5.41), but perhaps it would not be immediately obvious how to simplify it until you get (7.29) (or maybe it is: I did not try).
Long story short: for Riemannian submanifolds, both a direct computation as in (7.29) or using the general formula as in (5.41) work. The former requires you to know something about the connection; the latter requires you to know something about the Weingarten map. Of course, mathematically those are interchangeable (for submanifolds), but one may be more convenient than the other computationally.
About your second question, regarding submanifolds which are not Riemannian submanifolds: in those cases, we cannot use the simple general formulas as above, and we need to work from first principles: see the definitions of gradient and Hessian in Chs 3 and 5. Section 7.6 details how to do that for hyperbolic space (in the hyperboloid model, where the manifold in embedded in R^(n+1), and it does inherit a metric from R^(n+1) but it's the Minkowski metric, not a Euclidean one). It turns out that you can still pretty easily figure out the gradient, connection and Hessian in that case. For positive definite matrices, there are some results in Section 11.7, but few details for how to figure out the connection; there is a pointer to a paper by Sra & Hosseini, SIAM Opt. 2015, Section 3.
I hope this helps -- happy to answer more questions.