matlab - FDTD keeps overflowing on quasi boundaries and some other questions -

i have fdtd assignment, it's not going well. whenever put source in , let simulation run bit (note - if test have use debugger , step through loop increments otherwise there hefty overflow).

in of simulations i've seen, people break ex , ey separate parts, put them 1 large e matrix , used index offset (seen in ep(u,v+1) , e(u+1,v) calculate ey , ex independently).

the sources used reference were: http://fdtd.wikispaces.com/ , http://www.mathworks.com/matlabcentral/fileexchange/21000-tiny-fdtd-v1-0 1 acoustics, works well.

    close all;       %% user modifiable parameters     mu0 = pi*4e-7; % ph/µm     e0 =   8.854187e-12; % picofarads/micron     c = 1/sqrt(mu0*e0);      cellsizex = 100;      % size of yee-cell in microns     cellsizey = 100;      % size of yee-cell in microns     numx = 100; % number of cells in x direction     numy = 100; % number of cells in y direction     lambda = 700*10^-9;     dx = lambda/20;     dy = lambda/20;     dt =  (c*sqrt(dx^-2+dy^-2))^-1;      t0 = 100;       %index time of gaussian pulse peak     width = 10;     %peakyness of gaussian pulse      %% initialise h , e array     h = zeros(2*numx, 2*numy);            hp = zeros(2*numx, 2*numy);     e = zeros(2*numx+1,2*numy+1);         ep = zeros(2*numx+1,2*numy+1);         etemp = zeros(2*numx+1,2*numy+1);     htemp = zeros(2*numx, 2*numy);     p = zeros(2*numx+1,2*numy+1);     pp = zeros(2*numx+1,2*numy+1);      % scaling factors h , e fields      cex = dt/(dx*mu0);     cey = dt/(dy*mu0);     chx = dt/(dy*e0);     chy = dt/(dx*e0);     x = 2:2:2*numx; %2*numx-2;      %initialize permibilities     perm = ones(2*numx+1,2*numy+1);       %% fdtd loop n = 1:1500;      if(n < 200)         e(numx-10:2:numx+10,numy+1) =  1e-19*sin(2*pi*428e12*n*dt);    %exp(-.5 * ((n-t0)/width).^2); %insert hard source     end;      u = 2:2:2*numx-1         v = 2:2:2*numy-1             hp(u,v) = h(u,v) + (dt/mu0)*( -(e(u+1,v) - e(u-1,v))/dx) + ( e(u,v+1) - e(u,v-1)/dy ); % solving hplus             ep(u,v+1) = e(u,v+1) + (dt/(dy*e0))*(hp(u,v+2) - hp(u, v)); % solving ex plus             ep(u+1, v) = e(u+1, v) - (dt/(dx*e0))*(hp(u+2, v) - hp(u, v)); % solving ey plus          end     end;       % dirichlet boundary conditions     ep(1,:) = 0;     ep(:,1) = 0;     ep(2*numx+1,:) = 0;     ep(:,2*numx+1) = 0;      % plotting     surf(ep); shading interp; lighting phong; colormap hot; axis off; zlim([0 1]);     set(gcf,'color', [0 0 0], 'number', 'on', 'name', sprintf('fdtd project, step = %i', n));     %title(sprintf('time = %.2f microsec',n*dt*1e12),'color',[1 0 0],'fontsize', 22);      drawnow;      e = ep;     h = hp;    end;  end; 

2 things:

1 - ep(:,2*numx+1,:) = 0; correct? (i.e. should ep(:,2*numx+1) = 0;?)

2 - @ end, think may want change code to

if (max(max(ep)) > 1e-3)       e = ep/max(max(ep)); else e = ep; end if (max(max(hp)) > 1e-3)     h = hp/max(max(hp)); else h = hp; end 

as max(ep) give 1x100 matrix. 1e-3 protection, can lower whatever value choose.