# Finite Difference Algorithm on Non-Uniform Meshes for Modeling 2D Magnetotelluric Responses
## Appendix A MATLAB Program for Simulating TE-Mode Responses  
```
function[Ex,rho_a,phase]=MT2D_TE_FDM_Nonuniform(dy,dz_air,dz_earth,rho,fre) % Input arguments % dy: Mesh size in y direction % dz_air: Mesh size in z direction (air) % dz_earth: Mesh size in z direction (earth) % rho: Mesh resistivity % fre: Frequency % Output % Ex: Electric field % rho_a: Apparent resistivity % phase: Impedance phase mu=4e-7*pi; fre=logspace(-3,3,40); dz=[dz_air dz_earth]; Ny=length(dy); Nz_air=length(dz_air); Nz_earth=length(dz_earth); Nz=Nz_air+Nz_earth; L=sparse(Ny*Nz,Ny*Nz); R=sparse(Ny*Nz,1); for nf=1:1:size(fre,2) % Formed linear equation % Inner nodes for i=2:1:Nz-1  for j=2:1:Nz-1  for j=2:1:Ny-1  k=(j-1)*Nz+i;  L(k,k-Nz)=4/((dy(j)+dy(j+1))*(dy(j-1)+dy(j)));  L(k,k+Nz)=4/((dy(j)+dy(j+1))*(dy(j)+dy(j+1)));  L(k,k-1)=4/((dz(i)+dz(i+1))*(dz(i-1)+dz(i)));  L(k,k+1)=4/((dz(i)+dz(i+1))*(dz(i)+dz(i+1)));  L(k,k)=sqrt(-1)*z*pi*fre(nf)*mu/rho(i,j)-...  (4/((dy(j)+dy(j+1))*(dy(j-1)+dy(j)))+...  4/((dy(j)+dy(j+1))*(dy(j)+dy(j+1))))-...  (4/((dz(i)+dz(i+1))*(dz(i-1)+dz(i)))+...  4/((dz(i)+dz(i+1))*(dz(i)+dz(i+1)))); R(k,1)=0;  end end % Upper boundary i=1; for j=1:1:Ny  k=(j-1)*Nz+i;L(k,k)=l/dz(end)+sqrt(-sqrt(-1)*2*pi*fre(nf)*mu*(1/rho(i,j))); L(k,k-1)=-1/dz(end); R(k,l)=0; end % Left boundary for i=l:1:Nz for j=l:1:1:Ny k=(j-l)*Nz+i; if(j==l6&i>l6&i<Nz) L(k,k)=l;L(k,k+Nz)=-1; R(k,l)=0; end end % Right boundary for i=l:1:Nz for j=l:1:Ny k=(j-l)*Nz+i; if(j==Ny&&i>l6&i<Nz) L(k,k)=l;L(k,k-Nz)=-1; R(k,l)=0; end end end % Solving the linear equation u(:,nf)=L\R; u=full(u); u_new(:,:,nf)=reshape(u(:,nf),Nz,Ny); u1(:,nf)=u_new(Nz_air+1,:,nf); u2(:,nf)=u_new(Nz_air+2,:,nf); u3(:,nf)=u_new(Nz_air+3,:,nf); u4(:,nf)=u_new(Nz_air+4,:,nf); for i=l:1:Ny ux(i,nf)=(-l1*u1(i,nf)+18*u2(i,nf)-9*u3(i,nf)+2*u4(i,nf))... /(2*3*dz(Nz_air+1)); Zyx(i,nf)=u1(i,nf)/((1/(sqrt(-1)*2*pi*fre(nf)*mu))*ux(i,nf)); rho_a(i,nf)=abs(Zyx(i,nf))^2/(2*pi*fre(nf)*mu); phase(i,nf)=-atan(imag(Zyx(i,nf))/real(Zyx(i,nf)))*180/pi; end end Ex=u_new;

## Appendix B MATLAB Program for Simulating TM-Mode Responses

function [Hx,rho_a,phase]=MT2D_TM_FDM_Nonuniform(dy,dz,rho,fre) % Input arguments % dy: Mesh size in y direction % dz: Mesh size in z direction % rho: Mesh resistivity % fre: Frequency % Output arguments% Hx: Magnetic field % rho_a: Apparent resistivity % phase: Impedance phase mu=4e-7*pi; Ny=length(dy); Nz=length(dz); L=sparse(Ny*Nz,Ny*Nz); R=sparse(Ny*Nz,1); for nf=1:1:size(fre,2) % Formed linear equation % Inner nodes for i=2:1:Nz-1  for j=2:1:Ny-1  k=(j-1)*N2+i;  L(k,k-Nz)=(rho(i,j-1)+rho(i,j))*2/((dy(j)+dy(j+1))*(dy(j-1)+dy(j)));  L(k,k+Nz)=(rho(i,j+1)+rho(i,j))*2/((dy(j)+dy(j+1))*(dy(j)+dy(j+1)));  L(k,k-1)=(rho(i-1,j)+rho(i,j))*2/((dz(i)+dz(i+1))*(dz(i-1)+dz(i)));  L(k,k+1)=(rho(i+1,j)+rho(i,j))*2/((dz(i)+dz(i+1))*(dz(i)+dz(i+1)));  L(k,k)=sqrt(-1)*2*pi*free(nf)*mu-...  (rho(i,j-1)+rho(i,j))*2/((dy(j)+dy(j+1))*(dy(j-1)+dy(j)))...  (rho(i,j+1)+rho(i,j))*2/((dy(j)+dy(j+1))*(dy(j)+dy(j+1))))....  (rho(i-1,j)+rho(i,j))*2/((dz(i)+dz(i+1))*(dz(i-1)+dz(i)))...  (rho(i+1,j)+rho(i,j))*2/((dz(i)+dz(i+1))*(dz(i)+dz(i+1)));  R(k,1)=0;  end end % Upper boundary i=1; for j=1:1:Ny  k=(j-1)*Nz+i;  L(k,k)=1/dz(end)+sqrt(-sqrt(-1)*2*pi*fre(nf)*mu*(1/rho(i,j)));  L(k,k-1)=-1/dz(end);  R(k,1)=0; end % Left boundary for i=1:1:Nz  for j=1:1:Ny  k=(j-1)*N2+i;  if(j==1&&i>1&&i+Nz)  L(k,k)=1;L(k,k+Nz)=-1;  R(k,1)=0;  end end end
```end
* Right boundary for i=1:1:Nz for j=1:1:Ny k=(j-1)*Nz+i; if(j==Ny&i>1&5&i<Nz) L(k,k)=1;L(k,k-Nz)=-1; R(k,1)=0; end end end
* Solving the linear equation u(:,nf)=L\R; u=full(u); u_new(:,:,nf)=reshape(u(:,nf),Nz,Ny); u1(:,nf)=u_new(1,:,nf); u2(:,nf)=u_new(2,:,nf); u3(:,nf)=u_new(3,:,nf); u4(:,nf)=u_new(4,:,nf); for i=1:1:Ny u(i,nf)=(-11*u1(i,nf)+18*u2(i,nf)-9*u3(i,nf)+2*u4(i,nf))/(2+3*dz(1)); Zyx(i,nf)=rho(1,i)*ux(i,nf)/u1(i,nf); rho_a(i,nf)=abs(Zyx(i,nf))^2/(2*pi*fre(nf)*mu); phase(i,nf)=-atan(imag(Zyx(i,nf))/real(Zyx(i,nf)))*180/pi; end Hx=u_new;