This thesis presents the modeling of general medium constitutive relationships in the Transmission Line Matrix (TLM) method. The technique is shown for two- and three-dimensional cases. The procedure consists of decoupling the impulse scattering at the nodes from equations describing the medium. This is achieved by using nodal sources connected to the TLM node. The nodal sources are implemented with the state-variable description of the constitutive relationships. The technique requires only few modifications to the TLM algorithm. The procedure is validated for frequency-dependent, nonlinear, anisotropic and gyromagnetic media. This thesis also presents a dispersion analysis of TLM with frequency-dependent dielectrics. This study is performed in two- and three-dimensions by solving the dispersion relationship of TLM with nodal sources. The sources are used to model the frequency dependent dielectric. The study shows that the nodal source and stub-loaded models are equivalent for frequency independent dielectrics. The accuracy bounds of the TLM frequency-dependent dielectric model are presented. This thesis also investigates the physical origin of the coarseness and dispersion errors influencing two-dimensional TLM solutions of Maxwell's equations. The study is performed by solving the difference equations of the numerical method analytically. The results confirm a reduction of the accuracy of the discrete solution near field singularities. The solution of a partially filled waveguide is also investigated. The results show that TLM can have positive or negative frequency shifts, depending on the dielectric filling, excited mode and geometry. These results are also valid for the finite difference time domain method (FDTD).