I'm writing a program for finite difference time domain. Below, I wrote a program in Matlab to show FDTD E and H waves in 1 Dimension. What code do I have to add to make this a simulation in 2 Dimensions?

close all; clear all;

%define NN = 200; %# of iterations to perform N = 500; %# spatial samples in domain L = 5; %domain length in meters fs = 300e6; %source frequency in Hz ds = L/N; %spatial step in meters x = linspace(0,L,N); %x coordinate of spatial samples dt = ds/300e6; %"magic time step" eps0 = 8.854e-12; %permittivity of free space mu0 = pi*4e-7; %permeability of free space

%scale factors for E and H ae = ones(N,1)*dt/(ds*eps0); am = ones(N,1)*dt/(ds*mu0); as = ones(N,1); epsr = ones(N,1); mur= ones(N,1); sigma = zeros(N,1);

p = 1; for i=1:N epsr(i) = 1; mur(i) = 1; w1 = 0.5; w2 = 1.5; if (p==1) %dielectric window if (abs(x(i)-L/2)<0.5) epsr(i)=4; end end end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% ae = ae./epsr; am = am./mur; ae = ae./(1+dt*(sigma./epsr)/(2*eps0)); as = (1-dt*(sigma./epsr)/(2*eps0))./(1+dt*(sigma./epsr)/(2*eps0));

%initialize fields Hz = zeros(N,1); Ey = zeros(N,1);

for iter=1:NN %Yee Updataed Algorithm Ey(3) = Ey(3)+2*(1-exp(-((iter-1)/50)^2))*sin(2*pi*fs*dt*iter); %%%%%%%%%%%%%%%%% Hz(1) = Hz(2); %absorbing boundary conditions for H for i=2:N-1 %H field Update Hz(i) = Hz(i)-am(i)*(Ey(i+1)-Ey(i)); end %%%%%%%%%%%%%%% Ey(N) = Ey(N-1); %absorbing boundary conditions for E for i=2:N-1 % E field Update Ey(i) = as(i)*Ey(i)-ae(i)*(Hz(i)-Hz(i-1)); end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure(1) hold off plot(x,Ey,'b'); axis([3*ds L -2 2]); grid on; title('FDTD Gaussian Pulse'); hold on plot(x,380*Hz,'r');%multiplied by 380 so that the amplitude will reach same as E. legend('E Field','H-Field') xlabel('Time (Mintues)'); ylabel('Amplitude') pause(0); end