Finite Difference Time Domain (FDTD) - Python Help

Hey guys, I'm not sure if this is the correct place to post, if not I'm sorry!!

I'm trying to implement a 2D FDTD algorithm to solve Maxwell's Equations for a gaussian pulse source; but it's not working. The pulse explodes to ridiculous values and doesn't really propagate through the grid... I've gone through everything with a fine-toothed comb, but it's still broken :(

If there are any FDTD afficienardos out there who can scan my code to spot any obvious mistakes I'd be very grateful! Thanks :)

from __future__ import division import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D # SET UNITS kilohertz = 1e3 megahertz = 1e6 gigahertz = 1e9 # SET PHYSICAL CONSTANTS c0 = 299792458 eps0 = 8.85418782e-12 mu0 = 1.25663706e-6 # GRID PARAMETERS nmax = 1 # Air NLAM = 50 # Want to resolve the shortest wavelength (highest freq) with 50 points # SOURCE PARAMETERS fmax = 5.0 * gigahertz tau = 0.5 / fmax t0 = 5 * tau # SET COMPUTATIONAL CONSTANTS lam0 = c0 / fmax dx = lam0 / nmax / NLAM Nx = 200 dy = dx Ny = 200 steps = 1000 # CREATE GRID AXIS xa = np.linspace(0, Nx, Nx, endpoint = False) * dx ya = np.linspace(0, Ny, Ny, endpoint = False) * dy xa, ya = np.meshgrid(xa, ya) # INITIALISE MATERIALS TO FREE SPACE epsR = np.ones((Nx, Ny)) muR = np.ones((Nx, Ny)) # COMPUTE TIME STEP dt = min(dx, dy) / (2 * c0 ) # COMPUTE THE SOURCE time = np.linspace(0, steps, steps, endpoint = False) * dt source = np.exp(-((time - t0)/tau)**2) # COMPUTE UPDATE COEFFICIENTS mHx = - c0 * dt / muR mHy = c0 * dt / muR mDz = c0 * dt mEz = 1 / epsR # INITIALIZE FIELDS TO ZERO Hx = np.zeros((Nx, Ny)) Hy = np.zeros((Nx, Ny)) Dz = np.zeros((Nx, Ny)) Ez = np.zeros((Nx, Ny)) # MAIN FDTD LOOP for t in range(steps): # COMPUTE DZ for i in range(Nx): for j in range(Ny): if i == 0 and j == 0: Dz[0][0] = Dz[0][0] + mDz * ((Hy[0][0] - 0) / dx - (Hx[0][0] - 0) / dy) elif i == 0: Dz[0][j] = Dz[0][j] + mDz * ((Hy[0][j] - 0) / dx - (Hx[0][j] - Hx[0][j-1]) / dy) elif j == 0: Dz[i][0] = Dz[i][0] + mDz * ((Hy[i][0] - Hy[i-1][0]) / dx - (Hx[i][0] - 00) / dy) else: Dz[i][j] = Dz[i][j] + mDz * ((Hy[i][j] - Hy[i-1][j]) / dx - (Hx[i][j] - Hx[i][j-1]) / dy) # INJECT SOURCE Dz[int(round(Nx / 4))][int(round(Ny / 4))] += source[t] # COMPUTE EZ for i in range(Nx): for j in range(Ny): Ez[i][j] = mEz[i][j] * (Dz[i][j]) # COMPUTE HX for i in range(Nx): for j in range(Ny): if j == (Ny-1): Hx[i][Ny -1] = Hx[i][Ny -1] + mHx[i][Ny-1] * ((0 - Ez[i][Ny -1]) / dy) else: Hx[i][j] = Hx[i][j] + mHx[i][j] * ((Ez[i][j+1] - Ez[i][j]) / dy) # COMPUTE HY for i in range(Nx): for j in range(Ny): if i == (Nx-1): Hy[Nx-1][j] = Hy[Nx-1][j] + mHy[Nx-1][j] * (-(0 - Ez[Nx-1][j]) / dx) else: Hy[i][j] = Hy[i][j] + mHy[i][j] * (-(Ez[i+1][j] - Ez[i][j]) / dx) # CREATE FIGURE if t == 100: fig = plt.figure() ax = fig.add_subplot(111, projection='3d') ax.plot_wireframe(xa, ya, Ez, rstride = 5, cstride = 5) plt.savefig('C:\Users\Sam\Desktop\Python\FDTD\Images2D\Timestep %d' %t) plt.close()