The key to computation is visualization.

It appears that you are building a bunch of 'C???' arrays.  What are those?  Have you tried visualizing the data in those arrays after this section of code executes?  If that looks correct, the problem is elsewhere.  If they do not look correct, next visualize the arrays as they are being built.  This will give you a clue about what is going wrong.  You should probably look at one at a time instead of all at once.

Develop strong data visualization skills and you will find it a highly valuable tool for troubleshooting things like this.

Yes those are C arrays representing the constants multiplying the update equations in my code. The constants are delta t/(epsilon *delta x) which reduces to the courant number Sc, which is defined as 2-0.5 for 2D, multiplied by the impedance of free space imp0 divided by the relative permittivity or: (Sc*imp0)/EpsR.

They're defined in a header file for macros as:

#define ChyhG(G, M, N) G->chyh[(M) * SizeYG(G) + (N)]

#define ChyeG(G, M, N) G->chye[(M) * SizeYG(G) + (N)]

and my update equations are: 

void updateH2d(Grid *g){

        int mm, nn;

        if(Type == oneDGrid){

                for(mm = 0; mm < SizeX - 1; mm++)

                        Hy1(mm) = Chyh1(mm) * Hy1(mm) + Chye1(mm)*(Ez1(mm + 1) - Ez1(mm));

        } else {

                for (mm = 0; mm< SizeX; mm++)

                        for(nn = 0; nn < SizeY - 1;nn++)

                                Hx(mm, nn) = Chxh(mm, nn)* Hx(mm, nn) - Chxe(mm, nn) *(Ez(mm, nn + 1) - Ez(mm, nn));

                for(mm = 0; mm < SizeX - 1; mm++)

                        for(nn = 0; nn < SizeY; nn++)

                                Hy(mm, nn) = Chyh(mm, nn) * Hy(mm, nn) + Chye(mm, nn) * (Ez(mm + 1, nn) - Ez(mm, nn));

        }

        return;

}

/* update electric field */

void updateE2d(Grid *g) {

        int mm, nn;

        if (Type == oneDGrid){

                for (mm = 1; mm < SizeX - 1; mm ++)

                        Ez1(mm) = Ceze1(mm) * Ez1(mm) + Cezh1(mm) * (Hy1(mm) - Hy1(mm - 1));

        } else {

                for (mm = 1; mm < SizeX - 1; mm++)

                        for (nn = 1; nn < SizeY - 1; nn++)

                                Ez(mm, nn) = Ceze(mm, nn) * Ez(mm, nn) + Cezh(mm, nn) * ((Hy(mm, nn) - Hy(mm - 1, nn)) -(Hx(mm, nn) - Hx(mm, nn - 1)));

        }

        return;

}

I've plotted a grayscale image of the values for the constants as a function of their position and attached it below.

Also worth noting is my current problem trying to implement this advice is that there's some sort of infinite gain happening when I add in the layers. This doesn't happen when I take them out.

Are you incorporating loss?  It is easy to miss a sign and turn it into gain or to get something else wrong.  This may be what you are seeing so double check everything.  Even when simulating something with loss or gain, I highly recommend first running a simulation without that loss or gain and make sure your simulation obeys conservation of power.  That is, reflectance plus transmittance equals 100%.

R + T = 100%

Second, looking at your grid, I see a grating structure, but it is on the left side of your grid.  I am not sure how the rest of your code is working, but the grating may be interfering with how you inject a source.  That can also cause crazy things to happen.  Usually sources are incorporated into a part of the grid that is homogeneous.  Of course, this is not always the case but the source will need to be carefully considered when being injected where the materials are no homogeneous.

So the above image isn't a grating it's a quarter wave stack. The materials are dielectrics with refractive indexes of 1.46 and 2.1 respectively. The source comes in as a plane wave from a TFSF boundary on the left hand side of the image. The white in my image is all free space. So the source is injected on the left hand side of the image near the start of the white portion.

There is no gain or loss programed into my code. Any gain occurring is NOT intentional and something I'm trying to eliminate. 

This simulation could be done as a 1-D simulation (I have done so already). I'm now trying to get my results to agree in 2 and 3-D.

This is intended to be a simple simulation along the way towards testing my code for using it in future work. I haven't even added in conductivity yet. However adding two dielectrics has frustrated me for months and I've been unable to make any progress since May as a result. 

Have you tried simulating just air with your 2D code?

Have you tried setting one dimension to just 1 cell so that your 2D code effectively reduces to 1D so that you can compare data, update coefficients, etc?

Another thought...have you only incorporated your materials into some of the update coefficients?  For example, maybe your source is a linear polarization (say Ey), but you only incorporated your materials into the update coefficients for Ex.  If so, Ey never "feels" the materials.

I have simulated both vacuum on it's own as well as a single dielectric material in a vacuum successfully. I haven't tried the 1 cell with simulation in my 2D simulation yet. Right now I'm trying to put in periodic boundary conditions in the Z/Y directions to see if that changes anything (I have a question about boundary conditions I just asked!)

can you explain your latest post a little better. why would putting the material components in 3 dimensions effect what the other components "see"?

What I meant to say was that you have update coefficients for Ex, update coefficients for Ey, and update coefficients for Ez.  If your source is linearly polarized and only has the Ey component, for example, then the field would only see the materials that have been incorporated into the update coefficients for Ey.  If, for some reason, you incorporated your materials into the update coefficients for Ex and Ez and not Ey, then this could be the reason your fields are not seeing the materials.

I think I understand what is going on a little better. My question is this: what happens to a node at the boundary of two materials that is shared by two Yee cells that are in two different media. I've attached a picture to help explain, at nodes -1 and 6 the Ex nodes are in two different media, What permittivity should these nodes be set to and if they are set to one specific permittivity would that mean the Yee cell to the other side of the boundary now violates the divergence equations? If I want to code a structure that has a width of 35nm, as shown in the picture, how do I set the permittivities for the nodes as the boundaries? If the neighboring material also needs a specific width in the same dimension how would those node's permittivities then be set?

That is a great question.  The literature is full of "subcell" models to handle this.  An easy and good way to get through this is something I call the "2x grid technique."  Take a look at Lecture 12 here:

http://emlab.utep.edu/ee5390fdtd.htm

The short story is that you will calculate your grid, then calculate a grid that has twice as many points, then build your device on that 2x grid, then pull off the permittivity and permeability arrays from the 2x grid onto the 1x grid.  I hope the notes will make more sense out of what I said.  Once you understand what is really happening, you will be in a better position to learn more sophisticated ways of handling this.

Thanks! Could you point me to a place in the literature this is handled? I couldn't easily find it in Taflove or UFDTD? I'll also check out your lecture 

Lecture 12 here:

http://emlab.utep.edu/ee5390fdtd.htm

Thanks for the answer. So it's ok to have a Yee cell with different permittivities at different nodes without violating the divergence equation?

Yes.  For curved structures it is a necessity.  Also read through the blurring stuff.  That can help.