Perhaps all I need is to use the mean from predict and the covariance matrix from the undocumented function along with mvnrnd, similar to what's described in: https://www.mathworks.com/matlabcentral/answers/442492-sampling-based-on-a-correlation-matrix-for-multiple-parameters
Though I'm not sure if this is sampling from the posterior..
Thoughts?