From what I have seen, most people use Abaqus UEL, VUEL subroutines to do this if its not available in the software itself. So you only implement the element formulation, rest is taken care by the software.
Also many people use open source codes like dealii, mfem, feap, etc. to do this.
Thanks Mohammed abed Alruda and Akhilesh Pedgaonkar . I have created a software that formulates Shape functions for higher order Hexahedron (Lagrangian and serendipity), Tetrahedron and Custom element. I am sharing some screenshot of this software that is used to build a 22 node brick element.