3d FDTD code using MUR ABC

Please anyone help me.I can not solve 3d fdtd code using MUR ABC.I can not find any error in my program. So program output going to a large value.Herethis program.




format long
clear all
clc

%******************fundamenta constant*************************************


cc=299792458;
pi=3.14159;
epso=8.854e-12;
muo=4*pi*1.0e-7;
sig=0;
sigm=0;

%******************grid parameter***************************************** *
Nx=40;
Ny=40;
Nz=40;

Nxc=(Nx/2);
Nyc=(Ny/2);
Nzc=(Nz/2);

%*****************input value********************************************* *

freq=10e12;
ppw=15;
lambda=cc/freq;
%dx=lambda/ppw;
%dy=dx;
%dz=dy;
delta=.002;

dt=.002/(sqrt(3)*cc);
courant=(cc*dt)/.002;



%**********************Grid initilization**********************************

Ex=zeros(Nx-1,Ny,Nz);
cExe=zeros(Nx-1,Ny,Nz);
cExh=zeros(Nx-1,Ny,Nz);


Ey=zeros(Nx,Ny-1,Nz);
cEye=zeros(Nx,Ny-1,Nz);
cEyh=zeros(Nx,Ny-1,Nz);


Ez=zeros(Nx,Ny,Nz-1);
cEze=zeros(Nx,Ny,Nz-1);
cEzh=zeros(Nx,Ny,Nz-1);


Hx=zeros(Nx,Ny-1,Nz-1);
cHxh=zeros(Nx,Ny-1,Nz-1);
cHxe=zeros(Nx,Ny-1,Nz-1);


Hy=zeros(Nx-1,Ny,Nz-1);
cHyh=zeros(Nx-1,Ny,Nz-1);
cHye=zeros(Nx-1,Ny,Nz-1);


Hz=zeros(Nx-1,Ny-1,Nz);
cHzh=zeros(Nx-1,Ny-1,Nz);
cHze=zeros(Nx-1,Ny-1,Nz);


Eyoldleft=zeros(Nx,Ny-1,Nz);
Ezoldleft=zeros(Nx,Ny,Nz-1);
Eyoldright=zeros(Nx,Ny-1,Nz);
Ezoldright=zeros(Nx,Ny,Nz-1);
Exfrontold=zeros(Nx-1,Ny,Nz);
Ezfrontold=zeros(Nx,Ny,Nz-1);
Exrearold=zeros(Nx-1,Ny,Nz);
Ezrearold=zeros(Nx,Ny,Nz-1);
Exbottomold=zeros(Nx-1,Ny,Nz);
Eybottomold=zeros(Nx,Ny-1,Nz);
Extopold=zeros(Nx-1,Ny,Nz);
Eytopold=zeros(Nx,Ny-1,Nz);

%*********************electric field update coefficient********************

cExe(1:Nx-1,1:Ny,1:Nz)=(1-((sig*dt)/(2*epso)))/(1+((sig*dt)/(2*epso)));
cExh(1:Nx-1,1:Ny,1:Nz)=dt/((1+((sig*dt)/(2*epso)))*(epso*delta));
cEye(1:Nx,1:Ny-1,1:Nz)=(1-((sig*dt)/(2*epso)))/(1+((sig*dt)/(2*epso)));
cEyh(1:Nx,1:Ny-1,1:Nz)=dt/((1+((sig*dt)/(2*epso)))*(epso*delta));
cEze(1:Nx,1:Ny,1:Nz-1)=(1-((sig*dt)/(2*epso)))/(1+((sig*dt)/(2*epso)));
cEzh(1:Nx,1:Ny,1:Nz-1)=dt/((1+((sig*dt)/(2*epso)))*(epso*delta));

%*********************magnetic field update coefficient********************

cHxh(1:Nx,1:Ny-1,1:Nz-1)=(1-((sigm*dt)/(2*muo)))/(1+((sig*dt)/(2*muo)));
cHxe(1:Nx,1:Ny-1,1:Nz-1)=dt/((1+((sigm*dt)/(2*muo)))*(muo*delta));
cHyh(1:Nx-1,1:Ny,1:Nz-1)=(1-((sigm*dt)/(2*muo)))/(1+((sig*dt)/(2*muo)));
cHye(1:Nx-1,1:Ny,1:Nz-1)=dt/((1+((sigm*dt)/(2*muo)))*(muo*delta));
cHzh(1:Nx-1,1:Ny-1,1:Nz)=(1-((sigm*dt)/(2*muo)))/(1+((sig*dt)/(2*muo)));
cHze(1:Nx-1,1:Ny-1,1:Nz)=dt/((1+((sigm*dt)/(2*muo)))*(muo*delta));

%******************source parameter****************************************
t0 = 20.0;
spread = 6.0;
T = 0;

rtau=50.0e-12;
tau=rtau/dt;
ndelay=3*tau;
srcconst=-dt*3.0e+11;

fid=fopen('c:\\temp\\FDTDMURABC.txt','w');

%**************maximum time step*******************************************
nmax=500;
for n=1:nmax

T=T+1;
%***************************************Ex field**************************

for i=1:Nx-1
for j=2:Ny-1
for k=2:Nz-1

Ex(i,j,k)=cExe(i,j,k).* Ex(i,j,k)+ cExh(i,j,k).*((Hz(i,j,k)-...
Hz(i,j-1,k))-(Hy(i,j,k)-Hy(i,j,k-1)));

end
end
end
%**************************************Ey field****************************

for i=2:Nx-1
for j=1:Ny-1
for k=2:Nz-1

Ey(i,j,k)=cEye(i,j,k).* Ey(i,j,k)+cEyh(i,j,k).*((Hx(i,j,k)-...
Hx(i,j,k-1))-(Hz(i,j,k)-Hz(i-1,j,k)));

end
end
end
%*************************************Ez field****************************

for i=2:Nx-1
for j=2:Ny-1
for k=1:Nz-1

Ez(i,j,k)=cEze(i,j,k).* Ez(i,j,k)+cEzh(i,j,k).*((Hy(i,j,k)-...
Hy(i-1,j,k))-(Hx(i,j,k)-Hx(i,j-1,k)));

end
end
end
%..............................MUR ABC.....................................
m1 = (cc*dt - .002)/(cc*dt + .002);

% Left boundary

for j=1:Ny-1
for k=1:Nz
Ey(1,j,k) = Eyoldleft(2,j,k) + m1*(Ey(2,j,k) - Ey(1,j,k)); %left - Ey;
Eyoldleft(2,j,k)=Ey(2,j,k);

end
end

for j=1:Ny
for k=1:Nz-1
Ez(1,j,k) = Ezoldleft(2,j,k) + m1*(Ez(2,j,k) - Ez(1,j,k)); %left - Ez;
Ezoldleft(2,j,k)=Ez(2,j,k);
end
end

% Right boundary

for j=1:Ny-1
for k=1:Nz
Ey(Nx-1,j,k)= Eyoldright(Nx-2,j,k) + m1*(Ey(Nx-2,j,k) - Ey(Nx-1,j,k)); %right - Ey;
Eyoldright(Nx-2,j,k)=Ey(Nx-2,j,k);
end
end

for j=1:Ny
for k=1:Nz-1

Ez(Nx-1,j,k)= Ezoldright(Nx-2,j,k) + m1*(Ez(Nx-2,j,k) - Ez(Nx-1,j,k)); %right - Ez;
Ezoldright(Nx-2,j,k)=Ez(Nx-2,j,k);
end
end

% Front boundary
for i=1:Nx-1

for k=1:Nz
Ex(i, 1,k) = Exfrontold(i,2,k) + m1*(Ex(i,2,k) - Ex(i,1,k)); %front - Ex;
Exfrontold(i,2,k)=Ex(i,2,k);
end

end

for i=1:Nx

for k=1:Nz-1
Ez(i,1,k) = Ezfrontold(i,2,k) + m1*(Ez(i,2,k) - Ez(i,1,k)); %front - Ez;
Ezfrontold(i,2,k)=Ez(i,2,k);
end
end

% Rear boundary
for i=1:Nx-1
for k=1:Nz
Ex(i,Ny-1,k)= Exrearold(i,Ny-2,k) + m1*(Ex(i,Ny-2,k) - Ex(i,Ny-1,k)); % rear - Ex;
Exrearold(i,Ny-2,k)=Ex(i,Ny-2,k);
end
end

for i=1:Nx

for k=1:Nz-1
Ez(i,Ny-1,k)= Ezrearold(i,Ny-2,k) + m1*(Ez(i,Ny-2,k) - Ez(i,Ny-1,k)); % rear - Ey;
Ezrearold(i,Ny-2,k)=Ez(i,Ny-2,k);
end

end


% Bottom boundary
for i=1:Nx-1
for j=1:Ny

Ex(i,j,1) = Exbottomold(i,j,2) + m1*(Ex(i,j,2) - Ex(i,j,1)); % bottom - Ex;
Exbottomold(i,j,2)=Ex(i,j,2);
end
end

for i=1:Nx
for j=1:Ny-1

Ey(i,j,1) = Eybottomold(i,j,2) + m1*(Ey(i,j,2) - Ey(i,j,1)); % bottom - Ey;
Eybottomold(i,j,2)=Ey(i,j,2);

end
end

% Top boundary
for i=1:Nx-1
for j=1:Ny

Ex(i,j,Nz-1)= Extopold(i,j,Nz-2) + m1*(Ex(i,j,Nz-2) - Ex(i,j,Nz-1)); % top - Ex;
Extopold(i,j,Nz-2)=Ex(i,j,Nz-2);

end
end
for i=1:Nx
for j=1:Ny-1

Ey(i,j,Nz-1)= Eytopold(i,j,Nz-2) + m1*(Ey(i,j,Nz-2) - Ey(i,j,Nz-1)); % top - Ex;
Eytopold(i,j,Nz-2)=Ey(i,j,Nz-2);

end
end


%*************************************source****** ************************
pulse=exp(-0.5*((t0-T)/spread)^2);
Ez(Nxc,Nyc,Nzc)=pulse;

%Ez(Nxc,Nyc,1:Nzc)= srcconst*(n-ndelay)*exp(-((n-ndelay)^2/tau^2));

%Ez(Nxc,Nyc,Nzc)=Ez(Nxc,Nyc,Nz)+...
%srcconst*(n-ndelay)*exp(-((n-ndelay)^2/tau^2));

%**************************************Hx field****************************

for i=1:Nx
for j=1:Ny-1
for k=1:Nz-1

Hx(i,j,k)=cHxh(i,j,k)* Hx(i,j,k)+ cHxe(i,j,k)*((Ey(i,j,k+1)-...
Ey(i,j,k))-(Ez(i,j+1,k)-Ez(i,j,k)));

end
end
end

%****************************************Hy field**************************

for i=1:Nx-1
for j=1:Ny
for k=1:Nz-1

Hy(i,j,k)=cHyh(i,j,k)* Hy(i,j,k)+cHye(i,j,k)*((Ez(i+1,j,k)-...
Ez(i,j,k))-(Ex(i,j,k+1)-Ex(i,j,k)));

end
end
end

%*****************************************Hz field*************************

for i=1:Nx-1
for j=1:Ny-1
for k=1:Nz

Hz(i,j,k)=cHzh(i,j,k)* Hz(i,j,k)+cHze(i,j,k)*((Ex(i,j+1,k)-...
Ex(i,j,k))-(Ey(i+1,j,k)-Ey(i,j,k)));

end
end
end

fprintf( '%f\n',Ez(Nxc+2,Nyc,Nzc));



end




Plz hep me whats problem here.PLZ