Actually, it would be "overkill" to use Comsol instead of TMM for a multilayer system, but you certainly should arrive at the same answer if your meshing is fine enough and you use the same dielectric constants/indices of refraction for both methods (of course, there are a lot more errors you can stumble over when using Comsol). It is hard to make a remote diagnosis of where the problem might be situated. Unfortunately I don't know the book you cite, but if they use Abeles TMM everything should be fine. Can you reveal more about the materials you use, how many layers, what dielectric constants or indices of refractions etc.?
Since you are simulating a 1D photonic crystal in a waveguide, how are you simulating this with TMM? Are you using a one-dimensional TMM? If so, you will need accurate values of the effective refractive index for each layer. This could be where things are going wrong. Are you using a 2D or 3D TMM? If so, you should probably call that technique the method of lines or rigorous coupled-wave analysis. This is really what is needed to resolve the complexity of your problem I think.
There is a general class of numerical methods called the method of lines where all independent variable are solved numerically and the last is solved analytically. In electromagnetics, we interpret this as solving the cross section of a device numerically (i.e. calculating the modes that exist in the cross section) and then propagate them analytically in the longitudinal direction. In this manner, layer thicknesses in the longitudinal direction can be anything from nanometers to kilometers. Even more, very efficient techniques exist to connect multiple layers in the longitudinal direction.
I think “method of lines” is the better title for what you are calling transfer matrix methods. Transfer matrices is just one way that a method of lines code can connect layers stacked in the longitudinal direction. There are also scattering matrices, R matrices, H matrices, etc. In electromagnetics, we talk about the transfer matrix method (TMM), rigorous coupled-wave analysis (RCWA), Fourier model method (FMM), transfer matrix method with a plane wave basis, eigen-mode expansion technique, and the method of lines (MoL). All of these belong to the method of lines class of numerical methods. Most of the time when we say MoL in electromagnetics, for whatever reason, we are actually talking about the special case of method of lines where the finite-difference method is used in the transverse direction to calculate the modes. I think because MoL was popularized in electromagnetics by Reinhold Pregla who used mostly finite-differences.
With this background, the method I think you are asking about is called the method of lines based on finite-differences. The primary author for this method is Reinhold Pregla. A literature search will turn up a lot a great papers that he authored. I also teach this method in my Computational Electromagnetics course here in Lecture 24:
http://emlab.utep.edu/ee5390cem.htm
I teach this technique later in the course when the class has a lot of background already so the lecture is brief and probably not very understandable on its own. I suggest working through Lectures 2-6, 10-12, and then 24.
I am not clear on what it is you are doing. Are you calculating modes in a waveguide or are you simulating a device with a source? MoL can do everything, but it may be easier for you to use a different method if you are just calculating modes within a structure. For example, you could just use the simple finite-difference method described in Lecture 12 in my course. MoL is great for simulating devices, especially ones that are long in the longitudinal direction.
At this time, I want to caculate transmission, Reflection and absorption coefficient of Plasmonic waveguide with periodic structure. The structure consists of one Grahene sheet sandwiches between two SiO2 layers. Graphene doping, altered periodically. So we have 1D plasmonic crystal. Before, I simulated this structure by Comsol, I used a source current connected to Graphene sheet to excite plasmon. The extinction ratio is about 50 dB for 100 nm length. Then I used 1D TMM by using effective refractive index based on Bludov article "A PRIMER ON SURFACE PLADMON POLARITONS IN GRAPHENE" or Mier plasmonic book. TMM of Peter Markos(book) is very good. But 1D TMM results are very far from the Comsol Results, I am sure from Comsol. So, I should use MoL Based on your excellent comments to explain the Comsol. Also, 6 months ago, I simulate this structure by FDTD, I implemented a C code to simulate this structure. It is working good. I used a plane magnetic wave to excite Plasmon. But I don't know how to extract transmision and other coefficient from FDTD results!
I try to exactly implement your MoL lecture 24, for SPP Graphene 1D Photonic Crystal sandwiches between SiO2 Layers by MATLAB. I select Ky=0, and z is the direction of SPP wave(1D Photonic Crystal). First, Because of inhomogenity in x direction(Three Layers), epsilon gradian is not zero in two points. So, Could I use divergence formula to calculate Ez. May be, I should change the Ez calculation.
Second, I use dirichlet Boundary condition such as your thesis, page 97, Is it true?
Third, Why should we go to fourier domain from spatial domain, to calculate Ez.
Can I calculate the mean of Poynting vector in x direction in spatial domain to achieve reflection and transmission coefficients for a specific frequency.
Fourth, To Excite Plasmon, I used effective refractive index of Graphene waveguide in z direction for incident K=(neff)w/c, is it true? Also, I used Kx equality in Boundary Conditions and using dispersion relation in x direction to compute it, at all points.
I am not completely sure I understand the device you are simulating, but I think it is a 1D stack that supports an SPP and the direction of propagation is z. If you follow the formulation that is in the notes, you should not have to worry about doing anything special to handle a gradient in the permittivity. If you want to calculate Ez, finish the simulation and then use the formulas on Slide 10 in the latest version of the notes.
http://emlab.utep.edu/ee5390cem.htm
You can use Dirichlet boundary conditions as long as your SPP is not leaky. If it is, you must use a PML in addition to Dirichlet in order to absorb the outgoing waves.
The formulation presented is mostly for simulating diffraction through periodic structures. I chose to do this by calculating the power in each diffracted mode and then summing these numbers to get the overall reflectance and transmittance. The Fourier domain stuff is just an easy way to do this. For an SPP, it does not make sense to do it this way. You can integrate the z-component of the Poynting vector or get transmission and reflection or you can get this information from the global scattering matrix. You will need to identify the scattering coefficients that map to your SPP modes.
For MoL, you will have to calculate the eigen-modes in each cross section anyway. If I understand your problem, your source can just be one of these modes, the one that is the SPP. There is no need to actually calculate anything else and your source will be rigorous.