Generation, amplification, and evolution of magnetic fields in cool stars can be investigated by measuring the Zeeman effect in atomic and molecular lines observed in their spectra. In particular, Zeeman line broadening and polarization have been used for detecting magnetic fields in stellar atmospheres. Similar to the Sun, these fields are complex and height-dependent (i.e., comprise 3D structures) and require advanced diagnostics. Fortunately, many molecular lines dominating M-dwarf spectra, such as FeH, CaH, MgH, and TiO, are temperature- and Zeeman- sensitive and form at different atmospheric heights, which makes them excellent probes of magnetic fields on M dwarfs. Our goal is to analyze the complexity of magnetic fields in M dwarfs. We investigate how magnetic fields vary with the stellar temperature (i.e., mass) and how "surface" inhomogeneities are distributed in height - the dimension that is usually neglected in stellar magnetic studies. This is achieved by including many atomic and molecular species in our study. We have determined effective temperatures of the photosphere and of magnetic features, magnetic field strengths and filling factors for nine M dwarfs (M1-M7). Our chi^2^ analysis is based on a comparison of observed and synthetic intensity and circular polarization profiles (Stokes I and V) of many magnetically sensitive atomic and molecular lines in ten wavelength regions. Stokes profiles were calculated by solving polarized radiative transfer equations under the local thermodynamic equilibrium using model atmospheres. We have found that properties of magnetic structures depend on the analyzed atomic or molecular species and their formation heights within the atmosphere. Two types of magnetic features similar to those on the Sun have been found: one is cooler (starspots), while the other one is hotter (network, small-scale magnetic features). The magnetic field strength in both starspots and network is within 3kG to 6kG, on average it is 5kG for the M1-M7 spectral class range. These fields occupy a large fraction of M dwarf atmospheres at all heights, up to 100%. The plasma beta is less than one throughout the entire M dwarf atmospheres, implying that they are highly magnetized stars. A combination of many molecular and atomic species and a simultaneous analysis of intensity and circular polarization spectra have allowed us to better decipher the complexity of magnetic fields on M dwarfs, including their dependence on the height within the atmosphere. This work provides an opportunity to investigate a larger sample of M dwarfs as well as L-type brown dwarfs.