diff --git a/2D_FDTD_Aaron_Diego.py b/2D_FDTD_Aaron_Diego.py
new file mode 100644
index 0000000..1759c86
--- /dev/null
+++ b/2D_FDTD_Aaron_Diego.py
@@ -0,0 +1,83 @@
+import numpy as np
+import matplotlib.pyplot as plt
+#import scipy.constants
+import math
+from mpl_toolkits import mplot3d
+import matplotlib.animation as animation
+
+#Parámetros iniciales
+Ly =10
+Lx=10
+T=25
+K=500
+M=100
+N=100
+dx = Lx/M
+dy = Ly/N
+dt = T/K
+
+#Hacemos el espacio descreto
+lin_x = np.linspace(0,      Lx,        M+1, endpoint=True)
+lin_y = np.linspace(0,      Ly,        N+1, endpoint=True)
+X, Y = np.meshgrid(lin_x, lin_y)
+
+#Condiciones iniciales
+spread = 0.5
+H0 = np.exp(-(((X-5)**2)+(Y-5)**2)/spread)
+Ex0 = 0.0*X
+Ey0 = 0.0*Y
+
+#Hacemos la figura inicial
+fig = plt.figure()
+ax = plt.axes(projection='3d')
+ax.plot_surface(X, Y, H0, rstride=1, cstride=1,
+                cmap='viridis', edgecolor='none')
+ax.set_title('surface');
+ax.set_xlabel('Eje x')
+ax.set_ylabel('Eje y')
+ax.set_zlabel('H')
+
+#Declaramos variables de los campos
+ExOld = Ex0
+ExNew = np.zeros((len(X), len(Y)))
+EyOld = Ey0
+EyNew = np.zeros((len(X), len(Y)))
+probeH    = np.zeros(((K, len(X), len(Y))))
+hOld = H0
+hNew = np.zeros((len(X), len(Y)))
+
+#Calculamos con cada iteración los datos nuevos
+for n in range(K):
+    #Actualizar los datos
+    ExNew[1:-1,1:-1] = ExOld[1:-1,1:-1] + (dt/dy)*(hOld[1:-1,2:] - hOld[1:-1,:-2])
+    EyNew[1:-1,1:-1] = EyOld[1:-1,1:-1] - (dt/dx)*(hOld[2:,1:-1] - hOld[:-2,1:-1])
+    hNew[1:-1,1:-1]= hOld[1:-1,1:-1] + (ExNew[1:-1,2:] - ExNew[1:-1,:-2]) - (EyNew[2:,1:-1] - EyNew[:-2,1:-1])
+
+    #Condiciones de Frontera (reflectante)
+    ExNew[:,0] = ExOld[:,0] + (dt/dy)*(hOld[:,1] - hOld[:,N-1])
+    EyNew[0,:] = EyOld[0,:] - (dt/dx)*(hOld[1,:] - hOld[M-1,:])
+    ExNew[:,N] = ExOld[:,N] + (dt/dy)*(hOld[:,1] - hOld[:,N-1])
+    EyNew[M,:] = EyOld[M,:] - (dt/dx)*(hOld[1,:] - hOld[M-1,:])
+    hNew[0,:]= hOld[0,:] + (ExNew[0,:] - ExNew[0,:]) - (EyNew[1,:] - EyNew[M-1,:])
+    hNew[:,0]= hOld[:,0] + (ExNew[:,1] - ExNew[:,N-1]) - (EyNew[:,0] - EyNew[:,0])
+    hNew[M,:]= hOld[M,:] + (ExNew[M,:] - ExNew[M,:]) - (EyNew[1,:] - EyNew[M-1,:])
+    hNew[:,N]= hOld[:,N] + (ExNew[:,N] - ExNew[:,N-1]) - (EyNew[:,N] - EyNew[:,N])
+
+    #Actualizar los datos para la proxima iteración
+    ExOld[:,:] = ExNew[:,:]
+    EyOld[:,:] = EyNew[:,:]
+    hOld[:,:] = hNew[:,:]
+    probeH[n,:,:] = hNew[:,:]
+
+    #Hacer el plot
+def animate(i):
+
+    ax.collections.clear()
+    ax.plot_surface(X, Y, probeH[i], rstride=1, cstride=1,
+                cmap='viridis', edgecolor='none')
+    return
+
+anim = animation.FuncAnimation(fig, animate,
+                               frames=K, interval=50, blit=False)
+
+plt.show()
diff --git a/ETheo.py b/ETheo.py
new file mode 100644
index 0000000..7b4cf0f
--- /dev/null
+++ b/ETheo.py
@@ -0,0 +1,77 @@
+
+import numpy as np
+import math
+import scipy.constants
+import time
+import matplotlib.pyplot as plt
+import matplotlib.animation as animation
+
+c0   = scipy.constants.speed_of_light
+mu0  = scipy.constants.mu_0
+eps0 = scipy.constants.epsilon_0
+imp0 = math.sqrt(mu0 / eps0)
+
+L = 10
+dx        = 0.01
+finalTime = L/c0*2
+cfl       = .99
+
+dt = cfl * dx / c0
+numberOfTimeSteps = int( finalTime / dt )
+gridSize = int(L / dx)
+
+probeTime = np.zeros(numberOfTimeSteps) 
+probeE = np.zeros((gridSize,numberOfTimeSteps))
+gridE = np.zeros(gridSize)
+
+def F_analitica(x,t):
+    F_left = 0.5 * np.exp(-(x+1e9*t-L/2)**2)
+    F_right = 0.5 * np.exp(-(x-1e9*t-L/2)**2)
+    return F_left + F_right
+
+t = 0.0
+
+for n in range(numberOfTimeSteps):
+    probeTime[n] = t
+    x = 0.0
+    for i in range(gridSize):
+        probeE [i,n] = F_analitica (x,t)
+        gridE [i] = x
+        x += dx
+    t += dt
+
+
+# --- Creates animation ---
+fig = plt.figure(figsize=(8,4))
+ax1 = fig.add_subplot(1, 2, 1)
+ax1 = plt.axes(xlim=(gridE[0], gridE[-1]), ylim=(-1.1, 1.1))
+ax1.grid(color='gray', linestyle='--', linewidth=.2)
+ax1.set_xlabel('X coordinate [m]')
+ax1.set_ylabel('Field')
+line1,    = ax1.plot([], [], 'o', markersize=1)
+timeText1 = ax1.text(0.02, 0.95, '', transform=ax1.transAxes)
+
+# ax2 = fig.add_subplot(2, 2, 2)
+# ax2 = plt.axes(xlim=(gridE[0], gridE[-1]), ylim=(-1.1, 1.1))
+# ax2.grid(color='gray', linestyle='--', linewidth=.2)
+# ax2.set_xlabel('X coordinate [m]')
+# ax2.set_ylabel('Magnetic field [T]')
+# line2,    = ax2.plot([], [], 'o', markersize=1)
+# timeText2 = ax2.text(0.02, 0.95, '', transform=ax2.transAxes)
+
+def init():
+    line1.set_data([], [])
+    timeText1.set_text('')
+#    line2.set_data([], [])
+#    timeText2.set_text('')
+    return line1, timeText1#, line2, timeText2
+
+def animate(i):
+    line1.set_data(gridE, probeE[:,i])
+    timeText1.set_text('Time = %2.1f [ns]' % (probeTime[i]*1e9))
+#    line2.set_data(gridH, probeH[:,i]*100)
+#    timeText2.set_text('Time = %2.1f [ns]' % (probeTime[i]*1e9))
+    return line1, timeText1#, line2, timeText2
+
+anim = animation.FuncAnimation(fig, animate, init_func=init)
+plt.show()
\ No newline at end of file
diff --git a/FDTD_2D_480p.mov b/FDTD_2D_480p.mov
new file mode 100644
index 0000000..8ef3892
Binary files /dev/null and b/FDTD_2D_480p.mov differ
diff --git a/README.md b/README.md
new file mode 100644
index 0000000..7721ad7
--- /dev/null
+++ b/README.md
@@ -0,0 +1,4 @@
+# PyFDTD
+This program simulates the development of an electromagnetic wave with reflective boundary conditions based on Maxwell's equations.
+The electric component is the X-Y plane and the magnetic component is the Z-axis with time steps creating the animation.
+This was a project for my "Computational Methods in Theoretical Physics" class, part of the Physics and Mathematics master's program at the University of Granada