We present a coupled pair approach for studying few-body physics in harmonically trapped ultracold gases. The method is applied to a two-component Fermi system of $N$ particles. A stochastically variational gaussian expansion method is applied, focusing on optimization of the two-body correlations present in the strongly interacting, or unitary, limit. The groundstate energy of the four-, six- and eight-body problem with equal spin populations is calculated with high accuracy and minimal computational effort. We also calculate the structural properties of these systems and discuss their implication for the many-body ultracold gas and other few-body calculations.