Help with FDTD - need 2D FDTD C code

# include <math.h>
# include <stdlib.h>
# include <stdio.h>

# define IE 60
# define JE 60

int main ()
{
  float ga[IE][JE], dz[IE][JE], ez[IE][JE], hx[IE][JE], hy[IE][JE];
  int l, n, i, j, ic, jc, nsteps;
  float dx, dy, dt, T, epsz, pi, epsilon, sigma, mu, eaf;
  float t0, spread, pulse;
  FILE *fp,  *fopen();

  ic  =  IE/2;
  jc  =  JE/2;
  dx  =  .01;        /* cell size */
  dy  =  .01;        /* cell size */
  dt  =  dx/6e8;     /* Time steps */
  epsz  =  8.8e-12;
  pi  = 3.14159;
  mu  = 1.0e-7*4.0*pi;
  sigma  = 1.0e-8;

  for( j = 0; j < JE; j++) {
    printf("%2d", j);
    for(i = 0; i < IE; i++) {
      dz[i][j]  =  0.;
      ez[i][j]  =  0.;
      hx[i][j]  =  0.;
      hy[i][j]  =  0.;
      ga[i][j]  =  1.0;
      printf("%5.2f", ga[i][j]);
    }
    printf("\n");
  }

  t0 = 20.0;
  spread = 6.0;
  T = 0;
  nsteps = 1;

  while(nsteps > 0){
    printf("nsteps -> ");
    scanf("%d", &nsteps);
    printf("%d \n", nsteps);

    for(n = 1; n <= nsteps; n++) {
      T = T+1;

        /*  - Start of Main FDTD loop  - */

      /*  -  Put a Gaussian pulse in the middle  -  */
    
      pulse = sin(2*pi*1500*1e6*dt*T);
      ez[ic][jc] = pulse;

      /*  -  Calculate Ez field  -  */

     for(j = 1; j < JE; j++){
        for(i = 1; i < IE; i++) {
          ez[i][j] = (1.0 - sigma*dt/(2.0*epsz))/(1.0+sigma*dt)/(2.0*epsz)*ez[i][j]+(dt/epsz)
            /(1.0+sigma*dt/(2.0*epsz*(1/dx))*(hy[i][j] - hy[i - 1][j]) - (dt/epsz)/((1.0+sigma*dt)/(2.0*epsz))*(1/dy)
            *(hx[i][j] - hx[i][j - 1]));
        }
      }

      /*  -  Calculate Hx field  -  */
      for(j = 0; j < JE - 1; j++){
        for(i = 0; i < IE - 1; i++){
          hx[i][j] = hx[i][j] - (dt/mu*1/dy)*(ez[i][j+1] - ez[i][j]);
        }
      }
      
      /* Calculate Hy field */
      for(j = 0; j < JE - 1; j++){
        for(i = 0; i <= IE - 1; i++) {
          hy[i][j] = hy[i][j]+(dt/mu*1/dx)*(ez[i+1][j] - ez[i][j]);
        }
      }
    }

    /*  -  End of main FDTD loop  -  */

    for(j = 1; j < jc; j++) {
      printf("%2d", j);
      for(i = 1; i < ic; i++) {
        printf("%5.2f", ez[2*i][2*j]);
      }
      printf("\n");
    }
    printf("T = %5.0f \n", T);

    /* Write E field out to file "Ez" */

    fp = fopen("Ez", "w");
    for(j = 0; j < JE; j++) {
      for(i = 0; i < IE; i++) {
        fprintf(fp,"%d\t%d\t%6.3f",i, j, ez[i][j]);
      }
      fprintf(fp, "\n");
    }
    fclose(fp);
  }
}