Returning Jacobian in bvp4c/bvp5c

One can call bvp4c using something like this

options = bvpset('Stats', 'on', 'RelTol', 1e-4, 'AbsTol', 1e-8);
sol = bvp4c(mybvp, mybc, myinitial, options);

which gives a solution structure like

sol = 

          solver: 'bvp5c'
      parameters: 0.0501
               x: [1x259 double]
               y: [4x259 double]
           idata: [1x1 struct]
           stats: [1x1 struct]

I'd like to include the output of the Jacobian in the solution output (along with sol.x, sol.y)...like a sol.jac entry. Is there an easy way to do this? You're allowed to attach the Jacobian into the options structure, but I'd like the reverse to happen---i.e. a way to see the values that Matlab is using in it approximation of the Jaocbian.

[cwahl]

Has anyone been able to do this for bvp4c or bvp5c? If so, I would greatly appreciate if you took the time to post about how you managed to do that!