Unfortunately, helping you with BI-RME is not something I can do. I have not implemented that method, but I do think it is a great technique for what you are doing. I took a look at the papers and I understand your problem. You have a periodic waveguide which complicates the physics a little bit. If you are new to this sort of thing, I think it is going to take considerable effort to develop your own BI-RME algorithm. I suspect, however, it is only the answer to this problem that you seek and not necessarily developing BI-RME.
If you were to solve this with the finite-difference method, there are three approaches. Here they are in order of what I think would work best:
1. Finite-Difference Frequency-Domain: You would essentially build one 3D unit cell of your waveguide in memory and solve it as Eigen-value problem. This will give you what the field looks like as well as the propagation constant of your mode, which I think is really what you are ultimately after if I read the second paper correctly. While I only discuss this in 2D, you can see the basic method in Lectures 10-12 in my previous link.
2. Finite-Difference Time-Domain -- This is essentially the same technique, but in the time-domain. This will allow you to simulate a much larger device with great memory efficiency. I did not dig deep enough to see the size of your problem relative to wavelength, but since you are using metals it is probably rather small.
3. Method of Lines -- If you could make a simplifying assumption that your vertical rods have a square cross section, the MoL would be very fast and efficient for you. Otherwise, it will be slower to resolve the curved boundaries of the vertical rods.
Do you have access to any commercial software? I ask because the structure you are simulating is somewhat complex. Calculating the propagation constant of a standard waveguide is simple, but yours is periodic so it is a little more complicated.
I just thought of how I would do this...First I would model your top plate, bottom plate, and dielectric fill as a slab waveguide. This is a 1D simulation and it would give the effective refractive index of the mode. Second, I would run a 2D finite-difference simulation of a grid filled with the effective refractive index you just calculated and add the metal posts. This would give the overall solution. I think it would be a pretty accurate simulation as along as the top and bottom plates are close enough to have the substrate support only a single mode. I would try to reproduce the results from these papers to ensure it is working.
- - - Updated - - -
If you want to try this finite-difference approach, I would do the following:
1. Watch Lectures 10-12.
2. Look through the electronics notes for these same lectures. I have made some revisions, additions, and improvements.
3. Write a MATLAB code that simulates a slab waveguide. Here is a link about doing this you will find extremely helpful...Topic 7c Slab Waveguide Analysis
**broken link removed**
4. Now write a 2D waveguide calculator like what you saw in Step 2. Because you have a periodic waveguide, there is a slight twist to how this simulation will be done.