C program compiles fine but the ".exe" file stopped working, I need help in debugging the code

I have this c code which I have taken from a textbook, it compiles fine but when I run the ".exe" file it stops with the "stopped working" dialogue box. On debugging it I get the error "Unhandled exception at 0x01362289 in fdtd_in_free_space.exe: 0xC00000FD: Stack overflow (parameters: 0x00000000, 0x002C2000). occurred". I have windows 7 and I am using the 'Developer Command Prompt VS 2017' to compile and run the program. I can't figure out what is wrong in the code, I am new to programming so please forgive me if it is too simple.

I have looked at other similar questions here but the solutions did not help me.

/************************************
3D FDTD : DIPOLE IN FREE SPACE
Unhandled exception at 0x01362289 in fdtd_in_free_space.exe: 0xC00000FD: Stack overflow (parameters: 0x00000000, 0x002C2000). occurred
************************************/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#define XMAX 30
#define YMAX 30
#define ZMAX 30

int main() {

    float ga_x[XMAX][YMAX][ZMAX], ga_y[XMAX][YMAX][ZMAX], ga_z[XMAX][YMAX][ZMAX];
    float d_x[XMAX][YMAX][ZMAX], d_y[XMAX][YMAX][ZMAX], d_z[XMAX][YMAX][ZMAX];
    float e_x[XMAX][YMAX][ZMAX], e_y[XMAX][YMAX][ZMAX], e_z[XMAX][YMAX][ZMAX];
    float h_x[XMAX][YMAX][ZMAX], h_y[XMAX][YMAX][ZMAX], h_z[XMAX][YMAX][ZMAX];
    int i, j, k, l, n, x_mid, y_mid, z_mid, n_steps, n_pml;
    float dx, dt, T, eps_z, pi, eps, sigma, eaf;
    float t0, spread, pulse;
    FILE *fp, *fopen();

    x_mid = XMAX/2;
    y_mid = YMAX/2;
    z_mid = ZMAX/2;
    dx = 0.01;                      //cell size
    dt = dx/6e8;                    //time steps
    eps_z = 8.8e-12;
    pi = 3.14159;

    //Initializing the arrays
    for (k=0; k<ZMAX; k++){
        for (j=0; j<YMAX; j++){
            for (i=0; i<XMAX; i++){
                e_x[i][j][k] = 0.0;
                e_x[i][j][k] = 0.0;
                e_x[i][j][k] = 0.0;
                d_x[i][j][k] = 0.0;
                d_y[i][j][k] = 0.0;
                d_z[i][j][k] = 0.0;
                h_x[i][j][k] = 0.0;
                h_y[i][j][k] = 0.0;
                h_z[i][j][k] = 0.0;
                ga_x[i][j][k] = 0.0;
                ga_y[i][j][k] = 0.0;
                ga_z[i][j][k] = 0.0;
            }
        }
    }

    //specifying the dipole
    for (k=11; k<20; k++){
        ga_z[x_mid][y_mid][k] = 0.0;
    }
    ga_z[x_mid][y_mid][z_mid] = 0.0;

    //specifying the source pulse
    t0 = 20.0;
    spread = 6.0;

    T = 0;
    n_steps = 1;

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

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

            // calculate the D_x field
            for (k=1; k<ZMAX; k++){
                for (j=1; j<YMAX; j++){
                    for (i=1; i<XMAX; i++){
                        d_x[i][j][k] = d_x[i][j][k] + 0.5*(h_z[i][j][k] - h_z[i][j-1][k] - h_y[i][j][k] + h_y[i][j][k-1]);
                    }
                }
            }

            // calculate the D_y field
            for (k=1; k<ZMAX; k++){
                for (j=1; j<YMAX; j++){
                    for (i=1; i<XMAX; i++){
                        d_y[i][j][k] = d_y[i][j][k] + 0.5*(h_x[i][j][k] - h_x[i][j][k-1] - h_z[i][j][k] + h_y[i-1][j][k]);
                    }
                }
            }

            // calculate the D_z field
            for (k=1; k<ZMAX; k++){
                for (j=1; j<YMAX; j++){
                    for (i=1; i<XMAX; i++){
                        d_z[i][j][k] = d_z[i][j][k] + 0.5*(h_y[i][j][k] - h_y[i-1][j][k] - h_x[i][j][k] + h_x[i][j-1][k]);
                    }
                }
            }

            // Source
            pulse = exp(-0.5*(pow((t0-T)/spread,2.0)));
            d_z[x_mid][y_mid][z_mid] = pulse;

            // calculate the E from D field
            for (k=1; k<ZMAX-1; k++){
                for (j=1; j<YMAX-1; j++){
                    for (i=1; i<XMAX-1; i++){
                        e_x[i][j][k] = ga_x[i][j][k]*d_x[i][j][k];
                        e_y[i][j][k] = ga_y[i][j][k]*d_y[i][j][k];
                        e_z[i][j][k] = ga_z[i][j][k]*d_z[i][j][k];
                    }
                }
            }

            for (k=1; k<ZMAX-1; k++){
                for (j=1; j<YMAX-1; j++){
                    for (i=1; i<XMAX; i++){
                        h_x[i][j][k] = h_x[i][j][k] + 0.5*(e_y[i][j][k+1] - e_y[i][j][k] - e_z[i][j+1][k] + e_y[i][j][k]);
                    }
                }
            }

            for (k=1; k<ZMAX-1; k++){
                for (j=1; j<YMAX; j++){
                    for (i=1; i<XMAX-1; i++){
                        h_y[i][j][k] = h_y[i][j][k] + 0.5*(e_z[i+1][j][k] - e_z[i][j][k] - e_x[i][j][k+1] + e_x[i][j][k]);
                    }
                }
            }

            for (k=1; k<ZMAX; k++){
                for (j=1; j<YMAX-1; j++){
                    for (i=1; i<XMAX-1; i++){
                        h_z[i][j][k] = h_z[i][j][k] + 0.5*(e_x[i][j+1][k] - e_x[i][j][k] - e_y[i+1][j][k] + e_y[i][j][k]);
                    }
                }
            }
        }
    }

    // write the e-field out to a file
    fp = fopen("Ez", "w");
    for (j=0; j<YMAX; j++){
        for (i=0; i<XMAX; i++){
            fprintf(fp,"%7.4f",e_z[i][j][z_mid]);
        }
    }
    fclose(fp);

    return 0;
}