Yes, you are correct that the Jacobian calculated in parabolic/hyperbolic is only approximate. However, in my experience, ode15s performs very robustly with this approximate Jacobian. In general, ode15s prefers to cut the step size rather than request a new Jacobian. If the problem is well-posed, I have not encountered a case where ode15s was not able to complete a solution in hyperbolic and parabolic. And, of course, the Jacobian is simply used as a predictor-- it doesn't affect the actual solution.
Yes, pdenonlin does have some additional methods for calculating approximate Jacobians but they are all slower than the approach used in parabolic and hyperbolic.
In fact, there is some commented out code in parabolic that you can experiment with if you like. If you change the line
ode15s will use the numjac routine to calculate an approximate Jacobian numerically. In my experiments, this is much slower than the currently-implemented approach in parabolic. The solutions are the same, of course. One of the main strengths of ode15s is that it will terminate if it can't calculate a solution to the user-requested accuracy at any time step.