The life and death of massive stars (> 8 Msolar) have a significant role in the evolution of the universe through many cosmic phenomena. Massive stars are formed in the dense cores of the molecular clouds and the dense nature of the core (10^7 cm^-3) leads to a rich chemistry in these regions. These factors combined with their rarity and distance makes it difficult to probe the region to understand their evolution and unravel their complexity (Zinnecker & Yorke, 2007; Palla & Stahler, 1993). The use of complex organic molecules (COMs) to understand the evolution of MYSOs is an area of study that is quite new compared to other fields of astronomy and is the topic of this thesis. Currently, there has been some research into how COMs evolve during the evolution of these MYSOs but have mostly studied single sources with multiple COMs or less than ten sources with a few selected COMs. This thesis seeks to explore how COMs can be used to understand the evolution of MYSOs using a large sample of objects and whether definite evolutionary stages can be defined using COMs. The thesis also compares the COMs analysis with proposed formation pathways and three warm-up phase models for the evolution of COMs. This thesis details my PhD work which involved working with forty-one massive young stellar objects with comparable luminosities and distances. These sources have been observed with ALMA in Band 6. These observations covered the frequency range around 227GHz (in the LSB) and 241GHz (USB). The CASA software was used for the data reduction and data processing while the CASSIS software was used to fit the observations and to calculate the column density and excitation temperatures of the species towards the MYSOs. The process of data reduction using CASA and best fit parameter modelling using CASSIS are detailed in Chapter 2. Chapter 3 discusses the line analysis results obtained from the forty-one MYSOs. In this chapter, the characteristics and parameters of the COMs observed are discussed while Chapter 4 describes the results obtained for three of the MYSOs, which are representative of the other MYSOs observed. Chapter 5 compares the results obtained for the COMS observed in the MYSOs with other star-forming regions studied by researchers. Finally, Chapter 6 compares the COMs with each other using the Pearsonsâ€™ coefficient. The fractional abundances obtained for the COMs are also compared with the chemical warm-up models of star-formation. This study has clearly and distinctively found two stages of evolution of MYSOs for the forty-one MYSOs observed. An early, cold stage of evolution which is characterised by low temperatures (average temperature of 21 K), low abundances (average abundance of 6.54 x 10^-10 and a lack of gas-phase COMs except for CH3OH. The second stage of evolution is a later, warm stage of evolution which is characterised by high temperatures (average temperature of 231 K), high abundances (average abundance of 8.94 x 10^-7) and the presence of multiple COMs, this stage of evolution is triggered by the turn-on of the hot core phase of the MYSOs. There are a small number of MYSOs which showed a two-component temperature profile for CH3OH which could not be confidently placed into one of the two proposed evolutionary stages. These sources may indicate the existence of a third, short-lived evolutionary stage for MYSOs or they could alternatively be consistent with two distinct populations of gas along the line of sight to the MYSOs. This study can therefore confirm that COMs can be used to trace the evolution of MYSOs and also be used to show the distinct evolutionary stage of the MYSOs. The observation of COMs in MYSOs can also be used to confirm suggested formation pathways.