Ordinal[matin_]:=Module[{retmat,n,m,k},
n=Length[matin];
m=Length[matin[[1]]];
retmat=Table[0,{i,1,n},{j,1,m}];
k=1;
For[i=1,i<n+1,i=i+1,
For[j=1,j<m+1,j=j+1,
If[matin[[i,j]]==0,
retmat[[i,j]]=k;
k = k+1;
,
retmat[[i,j]]=0;
];
];
];
Return[retmat]
]

MBM[hdis_,bdis_,n_,m_]:=Module[{retmat,isubh,isubb,a,imin1,imax1,iminm,imaxm},
retmat=Table[0,{i,1,n},{j,1,m}];
isubh=Table[0,{j,1,m}];
isubb=Table[0,{j,1,m}];
For[j=1,j<m+1,j=j+1,isubh[[j]]=n-hdis[[j]]];
For[j=1,j<m+1,j=j+1,isubb[[j]]=n-bdis[[j]]];
For[j=1,j<m+1,j=j+1,For[i=1,i<n+1,i=i+1,If[isubh[[j]]≥i,retmat[[i,j]]=I,a=0]]];
For[j=1,j<m+1,j=j+1,For[i=1,i<n+1,i=i+1,If[isubb[[j]]≤i,retmat[[i,j]]=I,a=0]]];
imin1=isubh[[1]];
iminm=isubh[[m]];
imax1=isubb[[1]];
imaxm=isubb[[m]];
For[i=1,i<n+1,i=i+1,If[imin1<i<imax1,retmat[[i,1]]=I,a=0]];
For[i=1,i<n+1,i=i+1,If[iminm<i<imaxm,retmat[[i,m]]=I,a=0]];
Return[retmat]
]

Laplacian[hdis_,bdis_,n_,m_,psih_,dx_,dz_,boundarytype_]:=Module[{retmat,Dx,Dz,a,bk,dim,boundary,k,kplus,kminus,omat,bsubi,kb,psib},
If[boundarytype=="diffuse",a=0.0];
If[boundarytype=="reflective",a=2.0];
If[boundarytype=="absorbant",a=1.0];
boundary=MBM[hdis,bdis,n,m];
omat = Ordinal[boundary];
dim=Max[omat];
Dx=SparseArray[{i_,i_}->0.0,{dim,dim}];
Dz=SparseArray[{i_,i_}->0.0,{dim,dim}];
retmat=SparseArray[{i_,i_}->0.0,{dim,dim}];
bk=Table[0.0,{i,1,dim}];
psib=Table[0.0,{i,1,m}];
For[i=2,i<n,i=i+1,
For[j=2,j<m,j=j+1,
k=omat[[i,j]];
kplus=omat[[i,j+1]];
kminus=omat[[i,j-1]];
If[boundary[[i,j]]+boundary[[i,j-1]]+boundary[[i,j+1]]== 0,
Dx[[k,k]]=-2./dx^2;
Dx[[k,kplus]]=1./dx^2;
Dx[[k,kminus]]=1./dx^2;
];
k=omat[[i,j]];
kplus=omat[[i-1,j]];
kminus=omat[[i+1,j]];
If[boundary[[i,j]]+boundary[[i-1,j]]+boundary[[i+1,j]]==0,
Dz[[k,k]]=-2./dz^2;
Dz[[k,kplus]]=1./dz^2;
Dz[[k,kminus]]=1./dz^2;
];
If[boundary[[i-1,j]]==I && boundary[[i,j]]==0,
k=omat[[i,j]];
kminus=omat[[i+1,j]];
Dz[[k,k]]=-2./dz^2;
Dz[[k,kminus]]=1./dz^2;
bk[[k]]=-psih[[j]]
];
];
If[boundary[[i,2]]==0,
k=omat[[i,2]];
kplus=omat[[i,3]];
Dx[[k,k]]=-a/dx^2;
Dx[[k,kplus]]=a/dx^2;
];
If[boundary[[i,m-1]]==0,
k=omat[[i,m-1]];
kminus=omat[[i,m-2]];
Dx[[k,k]]=-a/dx^2;
Dx[[k,kminus]]=a/dx^2;
];
];
kb=omat[[n-bdis[[2]]-1,2]];
psib[[2]]=kb;
For[j=3,j<m,j=j+1,
If[bdis[[j]]==bdis[[j-1]],
kb=omat[[n-bdis[[j]]-1,j]];
psib[[j]]=kb;
];
If[bdis[[j]]==bdis[[j-1]]+1,
kb=omat[[n-bdis[[j]]-1,j-1]];
psib[[j]]=kb;
];
];
For[j=2,j<m,j=j+1,
k=omat[[n-bdis[[j]]-1,j]];
kplus=omat[[n-bdis[[j]]-2,j]];
Dz[[k,kplus]]=1./dz^2;
Dz[[k,k]]=-2./dz^2;
Dz[[k,kb]]=1./dz^2;
If[bdis[[j+1]]==bdis[[j]]+1,
kminus=omat[[n-bdis[[j]]-1,j-1]];
k=omat[[n-bdis[[j]]-1,j]];
Dx[[k,k]]=-2./dx^2;
Dx[[k,kminus]]=1./dx^2;
Dx[[k,kb]]=1./dx^2;
];
];
For[j=3,j<m-1,j=j+1,
If[hdis[[j-1]]==hdis[[j]]==hdis[[j+1]]+1,
k=omat[[n-hdis[[j]]+1,j]];
kminus=omat[[n-hdis[[j]]+1,j-1]];
Dx[[k,k]]=-2./dx^2;
Dx[[k,kminus]]=1./dx^2;
bk[[k]]=((-psih[[j+1]])/(dx^2))
];
If[hdis[[j-1]]+1==hdis[[j]]==hdis[[j+1]],
k=omat[[n-hdis[[j]]+1,j]];
kplus=omat[[n-hdis[[j]]+1,j+1]];
Dx[[k,k]]=-2./dx^2;
Dx[[k,kplus]]=1./dx^2;
bk[[k]]=((-psih[[j-1]])/(dx^2))
];
If[hdis[[j-1]]+1==hdis[[j]]==hdis[[j+1]]+1,
k=omat[[n-hdis[[j]]+1,j]];
Dx[[k,k]]=-2./dx^2;
bk[[k]]=((-(psih[[j-1]]+psih[[j+1]]))/(dx^2))
];
];
Return[{Dx+Dz,bk,omat}] 
]

TU[eta0_,phih0_,b_,n_,m_,dq_,ds_,Ns_,mod_]:=Module[{phihlast,phihc,etalast,etac,hdis,bdis,retmat,phitot,bnd,omat,Dmat,fxy,phi,r,phidq,j,l},
etalast=Table[eta0[(j-1)dq],{j,1,m}];
etac=etalast;
bdis=Table[Round[((b[(j-1)dq])/dq)],{j,1,m}];
hdis=Round[((1+etalast)/dq)];
phihlast=Table[phih0[(j-1)dq],{j,1,m}];
phihc=phihlast;
retmat=Table[0.0,{l,1,Ns+1},{j,1,m}];
phitot=Table[0.0,{l,1,Ns+1},{i,1,n},{j,1,m}];
For[l=1,l≤Ns&&Max[hdis]≤n &&Min[hdis]≥2,l=l+1,{Dmat,bnd,omat}=Laplacian[hdis,bdis,n,m,phihlast,dq,dq,"diffuse"];
fxy=LinearSolve[Dmat,bnd];
phi=Table[fxy[[omat[[n+1-i,j]]]],{i,1,n},{j,1,m}];
r=Table[phi[[hdis[[j]],j]],{j,1,m}];
For[j=3,j<m-1,j=j+1,
If[hdis[[j-1]]==hdis[[j]]==hdis[[j+1]],
phidq=(((phihlast[[j+1]]-phihlast[[j-1]]))/(2dq))];
If[hdis[[j-1]]==hdis[[j]]-1==hdis[[j+1]],
(((phihlast[[j+1]]-phihlast[[j-1]]))/(2dq))];
If[hdis[[j-1]]==hdis[[j]]==hdis[[j+1]]+1,
phidq=(((phihlast[[j]]-phihlast[[j-1]]+phihlast[[j+1]]-r[[j]]))/(2dq))];
If[hdis[[j-1]]==hdis[[j]]==hdis[[j+1]]-1,
phidq=(((r[[j+1]]-phihlast[[j-1]]))/(2dq))];
If[hdis[[j-1]]==hdis[[j]]+1==hdis[[j+1]],
phidq=(((r[[j+1]]-r[[j-1]]))/(2dq))];
If[hdis[[j-1]]+1==hdis[[j]]==hdis[[j+1]],
phidq=(((r[[j]]-phihlast[[j-1]]+phihlast[[j+1]]-phihlast[[j]]))/(2dq))];
If[hdis[[j-1]]-1==hdis[[j]]==hdis[[j+1]],
phidq=(((phihlast[[j+1]]-r[[j-1]]))/(2dq))];
If[hdis[[j-1]]-1==hdis[[j]]==hdis[[j+1]]+1,
phidq=(((phihlast[[j]]-r[[j-1]]+phihlast[[j+1]]-r[[j]]))/(2dq))];
If[hdis[[j-1]]+1==hdis[[j]]==hdis[[j+1]]-1,
phidq=(((r[[j]]-phihlast[[j-1]]+r[[j+1]]-phihlast[[j]]))/(2dq))];
(*time updating start*)
etac[[j]]=.5(etalast[[j+1]]+etalast[[j-1]])+ds(((phihlast[[j]]-r[[j]])/dq)-phidq(((etalast[[j+1]]-etalast[[j-1]]))/2));
phihc[[j]]=.5(phihlast[[j+1]]+phihlast[[j-1]])-ds(.5(phidq^2+(((phihlast[[j]]-r[[j]])/dq))^2)+dq hdis[[j]]-1);
(*time updating end*)
]
j=2;
phihc[[j]]=.5(phihlast[[j+1]]+phihlast[[j-1]])-ds(.5((((phihlast[[j+1]]-phihlast[[j-1]])/(2dq)))^2+(((phihlast[[j]]-r[[j]])/dq))^2)+dq hdis[[j]]-1);
etac[[j]]=etac[[j+1]];
j=m-1;
phihc[[j]]=.5(phihlast[[j+1]]+phihlast[[j-1]])-ds(.5((((phihlast[[j+1]]-phihlast[[j-1]])/(2dq)))^2+(((phihlast[[j]]-r[[j]])/dq))^2)+dq hdis[[j]]-1);
etac[[j]]=etac[[j-1]];
etac[[m]]=etac[[m-1]];
etac[[1]]=etac[[2]];
phihc[[1]]=2phihc[[2]]-phihc[[3]];
phihc[[m]]=2phihc[[m-1]]-phihc[[m-2]];
hdis=Round[(((1+etac))/dq)];
etalast=etac;
phihlast=phihc;
retmat[[l]]=etalast;
phitot[[l]]=phi;
If[Mod[l,mod]==0,Print[l]]
];
Return[{retmat,phitot}]
]