interpolation/approximation method for accurate evaluation of higher order derivatives of shape functions

Hello,

I  want to solve some problem from solid mechanics by means of meshless methods. Do you know some interpolation/approximation method which is able to accurate evaluate 3rd and 4th derivations of the shape functions and is not difficult to implement? MLS or PIM have problem with accuracy in these derivatives, as is known.
Thank you very much for your advices.