I developed a three dimensional beam model on Fortran with which I was hoping to obtain the displacements on the free edge of a model fixed on one end and free on the other (it suffers a perpendicular uniform load). However, as I increase the number of elements of the beam, the solution diverges from the theoretical solution and from the one obtained on Ansys for the same type of element. I don't think shear deformation could be a problem as I'm considering it on the stiffness matrix so what could be a possible explanation for this matter?