diff --git a/bin/RMT3d/Create_VTK.m b/bin/RMT3d/Create_VTK.m
new file mode 100644
index 0000000..1098ac5
--- /dev/null
+++ b/bin/RMT3d/Create_VTK.m
@@ -0,0 +1,78 @@
+function []=Create_VTK(x,y,z,value,VTK_file_name,type,cells)
+fid=fopen(VTK_file_name,'w');
+
+%% header
+header{1}='# vtk DataFile Version 3.0';
+header{2}=['VTK_3D_Pseudo test'];
+header{3}='ASCII';
+if strcmp(type,'Line')|| strcmp(type,'line')|| strcmp(type,'LINE')
+    header{4}='DATASET POLYDATA';
+else
+    header{4}='DATASET UNSTRUCTURED_GRID';
+end
+header{5}=['POINTS ',num2str(size(x,1)),' double'];
+for i=1:5
+    fprintf(fid,'%s\n',header{i});
+end
+%% point positions
+% for i=1:size(x,1)
+%     fprintf(fid,'%f %f %f\n',x(i),y(i),z(i));
+% end
+disp('writing positions')
+fprintf(fid,'%6.2f %6.2f %6.2f\n',[x,y,z]');
+% dlmwrite(VTK_file_name,[x y z],'-append','delimiter','\t','precision','%.1f')
+
+%% cells or lines
+fid=fopen(VTK_file_name,'a');
+if strcmp(type,'Line')|| strcmp(type,'line')|| strcmp(type,'LINE')
+    header{6}=['LINES ','1',' ',num2str(size(x,1)+1)];
+else
+    header{6}=['CELLS ',num2str(size(cells,1)),' ',num2str(5*size(cells,1))];
+end
+fprintf(fid,'%s\n',header{6});
+
+% cells(:,2)=[0:1:size(x,1)-1];
+Line_data=[size(x,1) 0:1:size(x,1)-1];
+if strcmp(type,'Line')|| strcmp(type,'line')|| strcmp(type,'LINE')
+    for i=1:size(x,1)+1
+        fprintf(fid,'%d\t',Line_data(i));
+    end
+else
+    %     for i=1:size(x,1)
+    %         fprintf(fid,'%d %d\n',cells(i,1),cells(i,2));
+    %     end
+    disp('writting cells')
+    cells=[4*ones(size(cells,1),1),cells];
+    fprintf(fid,'%i %i %i %i %i\n',cells');
+    %  dlmwrite(VTK_file_name,cells,'-append','delimiter','\t','precision','%i')
+end
+
+
+%% cells types
+if strcmp(type,'Line')||strcmp(type,'line')||strcmp(type,'LINE')
+    return
+end
+fid=fopen(VTK_file_name,'a');
+header{7}=['CELL_TYPES ',num2str(size(cells,1))];
+fprintf(fid,'%s\n',header{7});
+cells=10*ones(size(cells,1),1);
+%     for i=1:size(x,1)
+disp('writting cells type')
+%         fprintf(fid,'%d\n',cells(i,1));
+fprintf(fid,'%i\n',cells');
+%         dlmwrite(VTK_file_name,cells,'-append','delimiter','\t')
+%     end
+
+%% values
+fid=fopen(VTK_file_name,'a');
+header{8}=['CELL_DATA ',num2str(size(cells,1))];
+header{9}='SCALARS value double 1';
+header{10}='LOOKUP_TABLE default';
+for i=8:10
+    fprintf(fid,'%s\n',header{i});
+end
+disp('writting cell data')
+fprintf(fid,'%i\n',value');
+%  dlmwrite(VTK_file_name,value,'-append','delimiter','\t')
+fclose all;
+disp('done')
diff --git a/bin/RMT3d/calculate_layers.m b/bin/RMT3d/calculate_layers.m
new file mode 100644
index 0000000..ea8f2b5
--- /dev/null
+++ b/bin/RMT3d/calculate_layers.m
@@ -0,0 +1,6 @@
+
+function t=calculate_layers(n,t1,z)
+%% calculate layer thicknesses given nlayers and top and bottom
+exp_sum = @(x) ( x.^(n-1)-1)./(x-1)-z/t1;
+f = fzero(exp_sum,[1.0001 2]);
+t = t1*f.^(0:(n-2));
\ No newline at end of file
diff --git a/bin/RMT3d/em3d_MT_main.mlx b/bin/RMT3d/em3d_MT_main.mlx
new file mode 100644
index 0000000..482485f
Binary files /dev/null and b/bin/RMT3d/em3d_MT_main.mlx differ
diff --git a/bin/RMT3d/get_Boundary1DModel.mlx b/bin/RMT3d/get_Boundary1DModel.mlx
new file mode 100644
index 0000000..d251870
Binary files /dev/null and b/bin/RMT3d/get_Boundary1DModel.mlx differ
diff --git a/bin/RMT3d/get_BoundaryFieldAniEpsMu.mlx b/bin/RMT3d/get_BoundaryFieldAniEpsMu.mlx
new file mode 100644
index 0000000..3d5c33f
Binary files /dev/null and b/bin/RMT3d/get_BoundaryFieldAniEpsMu.mlx differ
diff --git a/bin/RMT3d/get_K_v2.mlx b/bin/RMT3d/get_K_v2.mlx
new file mode 100644
index 0000000..8cd00fb
Binary files /dev/null and b/bin/RMT3d/get_K_v2.mlx differ
diff --git a/bin/RMT3d/get_M_v2.mlx b/bin/RMT3d/get_M_v2.mlx
new file mode 100644
index 0000000..006b5ac
Binary files /dev/null and b/bin/RMT3d/get_M_v2.mlx differ
diff --git a/bin/RMT3d/get_abcd_volume.mlx b/bin/RMT3d/get_abcd_volume.mlx
new file mode 100644
index 0000000..92592eb
Binary files /dev/null and b/bin/RMT3d/get_abcd_volume.mlx differ
diff --git a/bin/RMT3d/get_boundary.mlx b/bin/RMT3d/get_boundary.mlx
new file mode 100644
index 0000000..ecf8ec5
Binary files /dev/null and b/bin/RMT3d/get_boundary.mlx differ
diff --git a/bin/RMT3d/get_data_MT.mlx b/bin/RMT3d/get_data_MT.mlx
new file mode 100644
index 0000000..6f045ca
Binary files /dev/null and b/bin/RMT3d/get_data_MT.mlx differ
diff --git a/bin/RMT3d/get_edges.m b/bin/RMT3d/get_edges.m
new file mode 100644
index 0000000..4173d8d
--- /dev/null
+++ b/bin/RMT3d/get_edges.m
@@ -0,0 +1,74 @@
+function [elems2edges, edges2nodes]=get_edges(elems2nodes)
+%function: [element2edges, edge2nodes]=get_edges(elems2nodes)
+%requires: deleterepeatedrows
+%generates edges of (triangular) triangulation defined in elems2nodes
+%elems2nodes is matrix, whose rows contain numbers of its element nodes 
+%element2edges returns edges numbers of each triangular element
+%edge2nodes returns two node numbers of each edge
+%example in 2D: [element2edges, edge2nodes]=get_edges([1 2 3; 2 4 3])
+%example in 3D: [element2edges, edge2nodes]=get_edges([1 2 3 4; 1 2 3 5; 1 2 4 6])
+ % Liu minghong 20221026   
+%2D case
+if (size(elems2nodes,2)==3)
+    %extracts sets of edges 
+    edges1=elems2nodes(:,[2 3]);
+    edges2=elems2nodes(:,[3 1]);
+    edges3=elems2nodes(:,[1 2]);
+
+    %as sets of their nodes (vertices)
+    vertices=zeros(size(elems2nodes,1)*3,2);
+    vertices(1:3:end,:)=edges1;
+    vertices(2:3:end,:)=edges2;
+    vertices(3:3:end,:)=edges3;
+
+    %repeated sets of nodes (joint edges) are eliminated 
+    [edges2nodes,elems2edges]=deleterepeatedrows(vertices);
+    elems2edges=reshape(elems2edges,3,size(elems2nodes,1))';
+end
+
+%3D case
+if (size(elems2nodes,2)==4)
+    %extracts sets of edges 
+    edges1=elems2nodes(:,[1 2]);
+    edges2=elems2nodes(:,[1 3]);
+    edges3=elems2nodes(:,[1 4]);
+    edges4=elems2nodes(:,[2 3]);     
+%     edges5=elems2nodes(:,[3 4]);
+%     edges6=elems2nodes(:,[4 2]);
+       % Liu minghong 20221026 
+    edges5=elems2nodes(:,[4 2]);
+    edges6=elems2nodes(:,[3 4]);
+        
+    %as sets of their nodes (vertices)
+    vertices=zeros(size(elems2nodes,1)*6,2);
+    vertices(1:6:end,:)=edges1;
+    vertices(2:6:end,:)=edges2;
+    vertices(3:6:end,:)=edges3;
+    vertices(4:6:end,:)=edges4;
+    vertices(5:6:end,:)=edges5;
+    vertices(6:6:end,:)=edges6;
+    
+    %repeated sets of nodes (joint edges) are eliminated 
+    [edges2nodes,elems2edges]=deleterepeatedrows(vertices);
+    elems2edges=reshape(elems2edges,6,size(elems2nodes,1))';
+    % *** 20221026
+    edges2nodes=sort(edges2nodes,2);
+    
+end
+
+
+function [matrix,I]=deleterepeatedrows(matrix)
+%function: [element2edges, edge2nodes]=edge_numbering(elements)
+%requires: deleterepeatedrows
+%generates edges of (triangular) triangulation defined in elements
+%elements is matrix, whose rows contain numbers of its element nodes 
+%element2edges returns edges numbers of each triangular element
+%edge2nodes returns two node numbers of each edge
+%example: [element2edges, edge2nodes]=edge_numbering([1 2 3; 2 4 3])
+
+
+%fast and short way suggested by John D'Ericco working in both 2D and 3D
+matrixs=sort(matrix,2);
+[dummy,J,I] = unique(matrixs,'rows');
+%I=reshape(I,size(matrixs,2),size(I,1)/size(matrixs,2));
+matrix=matrix(J,:);
\ No newline at end of file