Then, simply repeat what you do in the displacement method but use force/moment excitation instead to get flexibility.
Simply put.
You apply a load in position dof x, Fx1 and monitor displacements x1 and y1, which gives you in matrix notation the compliances C11 = x1/Fx1 and C12 = y1/Fx1. You then apply the load Fy1, monitor displacements and get compliances C21 = C12 = x1/Fy1 and C22 = y1/Fy1 .
You can assemble compliance into a matrix [C], its matrix inverse is the system stiffness [Z] = [C]-1. However & this is a BIG however, please note that a compliance matrix inverse produce stiffness only when all free dofs are part of the matrix inverse, i.e. for a general point you need 6 dofs (XYZRXRYRZ) and 6 dof generalized excitation (FxFyFzMxMyMz),
On a more practical note, when using FE, simply input unitary generalized excitation and you get the flexibility matrix terms automatically.