I want to write my own FEA program, when I calculate a cubic element(C3D8R) stiffness matrix, I found it is different from the matrix got from abaqus. Then, I found it is the hourglass that make it is different(I cancelled the hourglass, they are the same).  So my question is that, how to add the hourglass control into the Stiffness matrix by matlab code, sure, I just know how to add it is ok, code is not important. I only know artificial stiffness is ok, because i only use standard.I use matlab to write my FEA program, I want to know how to add hourglass into my program to get the same stiffness matrix got by matlab as abaqus C3D8R(default) stiffness matrix. I have read the paper Flanagan and Belytschko (1981) ,but i fail to understand its fortune code.