After solving for node values using the pde toolbox I use pdegrad to determine the gradient at all element centers and then use pdeprtni to interpolate the gradient values to the nodes. Is there a more accurate method to do the same in the pde toolbox. Thanks.