Yes, the structure has multiple regions with different dielectric constant. As a start, a dielectric region follows an empty region axially. So the epsilon matrix (or tensor if you prefer) changes its values as the axial index increases.

I've tried to alter my code, multiplying by the 1/sqrt(er) the timestep and the wavelength, using the largest available er, but the dt (timestep) isn't the same as without dielectric (empty guide),opposed to the theroretic approach where it should be the same. That's because numerical errors on rounding values and/or computer accuracy. So, I'm searching ways to overcome this issue but without sacryfying the underlying physical approach.An easy way is to multiply, by a factor

Yes, double precision

I've found the problem, it's about finding the smallest space step that satisfies the criterion space_step

Yes, I'm meaning that kind of errors, I already know that there are always rounding errors but, in a way some can be eliminated. That what I'm searching for, or a solution to a more-or-less same problem

Such as?Can you please let me know these methods? Because I would like to make my code as general as possible

I understood, I prefer not to change my equations, but use a more general formalism. Thanks for answering btw!

If your waveguide is filled axially with different axial layers of dielectrics you can use an exact analytical solution. Sucha asolution can be found in the papaers concerning axial mode matching method applied to dielectric resonators.

I agree with you Adam, but what I would like to see is the reflections. My question was about discretizing space. I had to use the smallest wavelength in order to build a right mesh for the simulation