diff --git a/example1/Input.m b/example1/Input.m
new file mode 100644
index 0000000..d4ad260
--- /dev/null
+++ b/example1/Input.m
@@ -0,0 +1,50 @@
+%%  SET PARAMETER
+addpath(genpath('D:\em3d_MT'));
+
+mesh.name  = 'mesh_layeredmodel'; % Contains grid information for the model, output by GeoMesh
+[mesh.node2coord,mesh.elem2node,mesh.entity] = load_comsol(mesh.name);
+% load_mesh;
+
+%% Set the conductivity
+P.airID     = [4];              % Domain ID of air
+P.sig_air   = 1e-8;             % Conductivity of air S/m
+P.sigp      = 1./[2000];        % Conductivity of background (residual domain) 
+P.sig_anoID = [2 1];            % Domain ID of the abnormal body
+P.sigax     = 1./[20000 20000]; % Conductivity of the abnormal body  
+P.sigay     = 1./[20000 20000];
+P.sigaz     = 1./[20000 20000];
+
+P.sig = Set_material_3Ani(mesh.entity,P.sig_anoID,P.airID,P.sigax,P.sigay,P.sigaz,P.sig_air,P.sigp) ;
+
+% Set the Euler rotation Angle of the abnormal body 
+P.aS   = [0 0];     
+P.aD   = [45 45];
+P.aL   = [0 0];
+P.AniAngle = Set_material_3Ani(mesh.entity,P.sig_anoID,P.airID,P.aS,P.aD,P.aL,0,0) ;
+P.sig  = get_sig_ani(P.sig(:,1),P.sig(:,2),P.sig(:,3),P.AniAngle(:,1),P.AniAngle(:,2),P.AniAngle(:,3));
+
+%% Set magnetic permitivity
+P.mu0 = 4*pi*1e-7 ;  % V·s/ (A·m)
+P.mu_anoID   = [];  % Domain ID of the abnormal body
+P.mura       = [];  % magnetic permitivity of the abnormal body
+P.mur = Set_material_3Ani(mesh.entity,P.mu_anoID,P.airID,P.mura,P.mura,P.mura,1,1) ;
+P.mur = P.mur(:,1);
+P.mu  = P.mur.*P.mu0 ;
+
+%% Set dielectric permittivity
+P.epsr_air   = 1;    
+P.epsrp      = 8;       % dielectric permittivity of background (residual domain) 
+P.eps_anoID  = [2 1];   % Domain ID of the abnormal body
+P.epsax      = [8 8];   % dielectric permittivity of the abnormal body
+P.epsay      = [5 5];
+P.epsaz      = [6 6];
+P.epsr = Set_material_3Ani(mesh.entity,P.eps_anoID,P.airID,P.epsax,P.epsay,P.epsaz,P.epsr_air,P.epsrp) ;
+P.epsr = get_sig_ani(P.epsr(:,1),P.epsr(:,2),P.epsr(:,3),P.AniAngle(:,1),P.AniAngle(:,2),P.AniAngle(:,3));
+P.eps0 = 8.854187817620389*10^-12 ;
+P.eps  = P.epsr *P.eps0 ;
+
+%% Set frequency
+freq  = logspace(4,6,11);
+omega = 2*pi*freq;
+nfreq = size(freq,2);
+
diff --git a/example1/RMT1Dexmple_ani_Esp_aD45_ChenHuang.mlx b/example1/RMT1Dexmple_ani_Esp_aD45_ChenHuang.mlx
new file mode 100644
index 0000000..c49b8f6
Binary files /dev/null and b/example1/RMT1Dexmple_ani_Esp_aD45_ChenHuang.mlx differ
diff --git a/example1/plot_freq_rho_1dvs3d.mlx b/example1/plot_freq_rho_1dvs3d.mlx
new file mode 100644
index 0000000..f7bf127
Binary files /dev/null and b/example1/plot_freq_rho_1dvs3d.mlx differ
diff --git a/example1/run.mlx b/example1/run.mlx
new file mode 100644
index 0000000..093fee0
Binary files /dev/null and b/example1/run.mlx differ