Long-range order evolution of self-assembled phosphonic acid multilayers as a function of temperature is studied here for two molecules with different alkyl chain length. By using synchrotron conventional diffraction, distinct order configurations are retrieved on phosphonic acid multilayers and their thermodynamic behavior monitored by energy-dispersive diffraction. This later technique allows us to observe the system behavior near order-disorder temperatures, as well as to determine the most stable configurations in the range from room temperature up to 120 degrees C. Planar order is also addressed by wide-angle X-ray scattering (WAXS) transmission experiments. Order parameter phase diagrams are built based on the experimental results, showing the dominant configuration at each temperature. The multilayer molecular long-range order retrieved from the experiments is corroborated by first principles calculations based on the Density Functional Theory. The bulk configurations depicted in this work are produced by molecule-molecule interactions and allow for future comparisons with the behavior of ordered molecules in few-monolayers configurations, commonly used in organic devices, where the presence of surfaces and interfaces strongly affects the molecule packing.