diff --git a/ExtractParameters.m b/ExtractParameters.m
new file mode 100644
index 0000000..8865d54
--- /dev/null
+++ b/ExtractParameters.m
@@ -0,0 +1,271 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%% Extract general parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% lattice parameter
+
+for i=1:length(DB.textdata(:,1))
+  idx=strcmp(DB.textdata{i,1},'a');
+  if idx==1
+    idx_a=i-1;
+    % break % removing the break makes it slower but more compatible between Matlab and Octave
+  end
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Epsilon static
+
+for i=1:length(DB.textdata(:,1))
+  idx=strcmp(DB.textdata{i,1},'Epsi_stat');
+  if idx==1
+    idx_Epsi=i-1;
+    % break % removing the break makes it slower but more compatible between Matlab and Octave
+  end
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% EP from the Kane model
+
+for i=1:length(DB.textdata(:,1))
+  idx=strcmp(DB.textdata{i,1},'EP_Kane');
+  if idx==1
+    idx_EP_K = i-1;
+    % break % removing the break makes it slower but more compatible between Matlab and Octave
+  end
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% EP from the Luttinger model
+
+% for i=1:length(DB.textdata(:,1))
+%   idx=strcmp(DB.textdata{i,1},'EP_Luttinger');
+%   if idx==1
+%     idx_EP_L = i-1;
+%     % break % removing the break makes it slower but more compatible between Matlab and Octave
+%   end
+% end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Energy level of the band 7c
+
+%for i=1:length(DB.textdata(:,1))
+%  idx=strcmp(DB.textdata{i,1},'Eg7c');
+%  if idx==1
+%    idx_Eg7c = i-1;
+%    % break % removing the break makes it slower but more compatible between Matlab and Octave
+%  end
+%end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Energy level of the band 8c
+
+%for i=1:length(DB.textdata(:,1))
+%  idx=strcmp(DB.textdata{i,1},'Eg8c');
+%  if idx==1
+%    idx_Eg8c = i-1;
+%    % break % removing the break makes it slower but more compatible between Matlab and Octave
+%  end
+%end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Energy level of the band 6v
+
+%for i=1:length(DB.textdata(:,1))
+%  idx=strcmp(DB.textdata{i,1},'Eg6v');
+%  if idx==1
+%    idx_Eg6v = i-1;
+%    % break % removing the break makes it slower but more compatible between Matlab and Octave
+%  end
+%end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Energy level of the band 6c (the band gap in a direct band gap semiconductor)
+
+for i=1:length(DB.textdata(:,1))
+  idx=strcmp(DB.textdata{i,1},'Eg6cG');
+  if idx==1
+    idx_Eg6c = i-1;
+    % break % removing the break makes it slower but more compatible between Matlab and Octave
+  end
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Bandgap Temperature dependency parameter "alpha" at the Gamma point
+
+for i=1:length(DB.textdata(:,1))
+  idx=strcmp(DB.textdata{i,1},'alphaG');
+  if idx==1
+    idx_alphaG = i-1;
+    % break % removing the break makes it slower but more compatible between Matlab and Octave
+  end
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Bandgap Temperature dependency parameter "beta" at the Gamma point
+
+for i=1:length(DB.textdata(:,1))
+  idx=strcmp(DB.textdata{i,1},'betaG');
+  if idx==1
+    idx_betaG = i-1;
+    % break % removing the break makes it slower but more compatible between Matlab and Octave
+  end
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% split-off band energy Delta-so
+
+for i=1:length(DB.textdata(:,1))
+  idx=strcmp(DB.textdata{i,1},'Dso');
+  if idx==1
+    idx_Dso = i-1;
+    % break % removing the break makes it slower but more compatible between Matlab and Octave
+  end
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Luttinger parameter for the electrons
+
+% for i=1:length(DB.textdata(:,1))
+%   idx=strcmp(DB.textdata{i,1},'F');
+%   if idx==1
+%     idx_F = i-1;
+%     % break % removing the break makes it slower but more compatible between Matlab and Octave
+%   end
+% end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Luttinger parameter gamma1 
+
+% for i=1:length(DB.textdata(:,1))
+%   idx=strcmp(DB.textdata{i,1},'gamma1');
+%   if idx==1
+%     idx_g1 = i-1;
+%     % break % removing the break makes it slower but more compatible between Matlab and Octave
+%   end
+% end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Luttinger parameter gamma2
+
+% for i=1:length(DB.textdata(:,1))
+%   idx=strcmp(DB.textdata{i,1},'gamma2');
+%   if idx==1
+%     idx_g2 = i-1;
+%     % break % removing the break makes it slower but more compatible between Matlab and Octave
+%   end
+% end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Luttinger parameter gamma3
+
+% for i=1:length(DB.textdata(:,1))
+%   idx=strcmp(DB.textdata{i,1},'gamma3');
+%   if idx==1
+%     idx_g3 = i-1;
+%     % break % removing the break makes it slower but more compatible between Matlab and Octave
+%   end
+% end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Valence band offset
+
+for i=1:length(DB.textdata(:,1))
+  idx=strcmp(DB.textdata{i,1},'VBO');
+  if idx==1
+    idx_VBO = i-1;
+    % break % removing the break makes it slower but more compatible between Matlab and Octave
+  end
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% mass for the 1band Schrodinger solver
+
+for i=1:length(DB.textdata(:,1))
+  idx=strcmp(DB.textdata{i,1},'me');
+  if idx==1
+    idx_me = i-1;
+    % break % removing the break makes it slower but more compatible between Matlab and Octave
+  end
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%% Strain parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% ac Conduction band Strain parameter
+
+for i=1:length(DB.textdata(:,1))
+  idx=strcmp(DB.textdata{i,1},'ac');
+  if idx==1
+    idx_ac = i-1;
+    % break % removing the break makes it slower but more compatible between Matlab and Octave
+  end
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% av Valence band Strain parameter
+
+for i=1:length(DB.textdata(:,1))
+  idx=strcmp(DB.textdata{i,1},'av');
+  if idx==1
+    idx_av = i-1;
+    % break % removing the break makes it slower but more compatible between Matlab and Octave
+  end
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% bv Valence band Strain parameter
+
+for i=1:length(DB.textdata(:,1))
+  idx=strcmp(DB.textdata{i,1},'bv');
+  if idx==1
+    idx_bv = i-1;
+    % break % removing the break makes it slower but more compatible between Matlab and Octave
+  end
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% dv Valence band Strain parameter
+
+%for i=1:length(DB.textdata(:,1))
+%  idx=strcmp(DB.textdata{i,1},'dv');
+%  if idx==1
+%    idx_dv = i-1;
+%    % break % removing the break makes it slower but more compatible between Matlab and Octave
+%  end
+%end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% c11 Valence band Strain parameter
+
+for i=1:length(DB.textdata(:,1))
+  idx=strcmp(DB.textdata{i,1},'c11');
+  if idx==1
+    idx_c11 = i-1;
+    % break % removing the break makes it slower but more compatible between Matlab and Octave
+  end
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% c12 Valence band Strain parameter
+
+for i=1:length(DB.textdata(:,1))
+  idx=strcmp(DB.textdata{i,1},'c12');
+  if idx==1
+    idx_c12 = i-1;
+    % break % removing the break makes it slower but more compatible between Matlab and Octave
+  end
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% c44 Valence band Strain parameter
+
+%for i=1:length(DB.textdata(:,1))
+%  idx=strcmp(DB.textdata{i,1},'c44');
+%  if idx==1
+%    idx_c44 = i-1;
+%    % break % removing the break makes it slower but more compatible between Matlab and Octave
+%  end
+%end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\ No newline at end of file
diff --git a/Library.m b/Library.m
new file mode 100644
index 0000000..821101b
--- /dev/null
+++ b/Library.m
@@ -0,0 +1,709 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Library load %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+DB = importdata('materialDB_ZB.csv',',');
+%DB = importdata('materialDB_WZ.csv',',');
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%% Here is a patch to be able to load the tables in Matlab AND Octave %%%%%%
+% Matlab see the header in multiple cells while Octave see the header in one cell only
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+M='GaAs';
+if length(DB.textdata(1,:))==1  %% Octave data load
+
+    DB.textdata{1,1}=[DB.textdata{1,1} ',']; % patch, add a comma "," at the end
+    idxM=strfind(DB.textdata{1,1},',');
+    idx=strfind(DB.textdata{1,1},[',' M ',']);
+    idxM=find(idxM==idx);
+    
+    GaAs = DB.data(:,idxM)';
+else  %% Matlab data load
+
+    for i=1:length(DB.textdata(1,:))
+      idx=strcmp(DB.textdata{1,i},M);
+      if idx==1
+        GaAs=DB.data(:,i-1)';
+        % break % removing the break makes it slower but more compatible between Matlab and Octave
+      end
+    end 
+    
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+M='AlAs';
+if length(DB.textdata(1,:))==1  %% Octave data load
+
+    DB.textdata{1,1}=[DB.textdata{1,1} ',']; % patch, add a comma "," at the end
+    idxM=strfind(DB.textdata{1,1},',');
+    idx=strfind(DB.textdata{1,1},[',' M ',']);
+    idxM=find(idxM==idx);
+    
+    AlAs = DB.data(:,idxM)';
+else  %% Matlab data load
+
+    for i=1:length(DB.textdata(1,:))
+      idx=strcmp(DB.textdata{1,i},M);
+      if idx==1
+        AlAs=DB.data(:,i-1)';
+        % break % removing the break makes it slower but more compatible between Matlab and Octave
+      end
+    end 
+    
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+M='InAs';
+if length(DB.textdata(1,:))==1  %% Octave data load
+
+    DB.textdata{1,1}=[DB.textdata{1,1} ',']; % patch, add a comma "," at the end
+    idxM=strfind(DB.textdata{1,1},',');
+    idx=strfind(DB.textdata{1,1},[',' M ',']);
+    idxM=find(idxM==idx);
+    
+    InAs = DB.data(:,idxM)';
+else  %% Matlab data load
+
+    for i=1:length(DB.textdata(1,:))
+      idx=strcmp(DB.textdata{1,i},M);
+      if idx==1
+        InAs=DB.data(:,i-1)';
+        % break % removing the break makes it slower but more compatible between Matlab and Octave
+      end
+    end 
+    
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+M='InSb';
+if length(DB.textdata(1,:))==1  %% Octave data load
+
+    DB.textdata{1,1}=[DB.textdata{1,1} ',']; % patch, add a comma "," at the end
+    idxM=strfind(DB.textdata{1,1},',');
+    idx=strfind(DB.textdata{1,1},[',' M ',']);
+    idxM=find(idxM==idx);
+    
+    InSb = DB.data(:,idxM)';
+else  %% Matlab data load
+
+    for i=1:length(DB.textdata(1,:))
+      idx=strcmp(DB.textdata{1,i},M);
+      if idx==1
+        InSb=DB.data(:,i-1)';
+        % break % removing the break makes it slower but more compatible between Matlab and Octave
+      end
+    end 
+    
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+M='GaSb';
+if length(DB.textdata(1,:))==1  %% Octave data load
+
+    DB.textdata{1,1}=[DB.textdata{1,1} ',']; % patch, add a comma "," at the end
+    idxM=strfind(DB.textdata{1,1},',');
+    idx=strfind(DB.textdata{1,1},[',' M ',']);
+    idxM=find(idxM==idx);
+    
+    GaSb = DB.data(:,idxM)';
+else  %% Matlab data load
+
+    for i=1:length(DB.textdata(1,:))
+      idx=strcmp(DB.textdata{1,i},M);
+      if idx==1
+        GaSb=DB.data(:,i-1)';
+        % break % removing the break makes it slower but more compatible between Matlab and Octave
+      end
+    end 
+    
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+M='AlSb';
+if length(DB.textdata(1,:))==1  %% Octave data load
+
+    DB.textdata{1,1}=[DB.textdata{1,1} ',']; % patch, add a comma "," at the end
+    idxM=strfind(DB.textdata{1,1},',');
+    idx=strfind(DB.textdata{1,1},[',' M ',']);
+    idxM=find(idxM==idx);
+    
+    AlSb = DB.data(:,idxM)';
+else  %% Matlab data load
+
+    for i=1:length(DB.textdata(1,:))
+      idx=strcmp(DB.textdata{1,i},M);
+      if idx==1
+        AlSb=DB.data(:,i-1)';
+        % break % removing the break makes it slower but more compatible between Matlab and Octave
+      end
+    end 
+    
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+M='InP';
+if length(DB.textdata(1,:))==1  %% Octave data load
+
+    DB.textdata{1,1}=[DB.textdata{1,1} ',']; % patch, add a comma "," at the end
+    idxM=strfind(DB.textdata{1,1},',');
+    idx=strfind(DB.textdata{1,1},[',' M ',']);
+    idxM=find(idxM==idx);
+    
+    InP = DB.data(:,idxM)';
+else  %% Matlab data load
+
+    for i=1:length(DB.textdata(1,:))
+      idx=strcmp(DB.textdata{1,i},M);
+      if idx==1
+        InP=DB.data(:,i-1)';
+        % break % removing the break makes it slower but more compatible between Matlab and Octave
+      end
+    end 
+    
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+M='GaP';
+if length(DB.textdata(1,:))==1  %% Octave data load
+
+    DB.textdata{1,1}=[DB.textdata{1,1} ',']; % patch, add a comma "," at the end
+    idxM=strfind(DB.textdata{1,1},',');
+    idx=strfind(DB.textdata{1,1},[',' M ',']);
+    idxM=find(idxM==idx);
+    
+    GaP = DB.data(:,idxM)';
+else  %% Matlab data load
+
+    for i=1:length(DB.textdata(1,:))
+      idx=strcmp(DB.textdata{1,i},M);
+      if idx==1
+        GaP=DB.data(:,i-1)';
+        % break % removing the break makes it slower but more compatible between Matlab and Octave
+      end
+    end 
+    
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+M='AlP';
+if length(DB.textdata(1,:))==1  %% Octave data load
+
+    DB.textdata{1,1}=[DB.textdata{1,1} ',']; % patch, add a comma "," at the end
+    idxM=strfind(DB.textdata{1,1},',');
+    idx=strfind(DB.textdata{1,1},[',' M ',']);
+    idxM=find(idxM==idx);
+    
+    AlP = DB.data(:,idxM)';
+else  %% Matlab data load
+
+    for i=1:length(DB.textdata(1,:))
+      idx=strcmp(DB.textdata{1,i},M);
+      if idx==1
+        AlP=DB.data(:,i-1)';
+        % break % removing the break makes it slower but more compatible between Matlab and Octave
+      end
+    end 
+    
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+M='Si';
+if length(DB.textdata(1,:))==1  %% Octave data load
+
+    DB.textdata{1,1}=[DB.textdata{1,1} ',']; % patch, add a comma "," at the end
+    idxM=strfind(DB.textdata{1,1},',');
+    idx=strfind(DB.textdata{1,1},[',' M ',']);
+    idxM=find(idxM==idx);
+    
+    Si = DB.data(:,idxM)';
+else  %% Matlab data load
+
+    for i=1:length(DB.textdata(1,:))
+      idx=strcmp(DB.textdata{1,i},M);
+      if idx==1
+        Si=DB.data(:,i-1)';
+        % break % removing the break makes it slower but more compatible between Matlab and Octave
+      end
+    end 
+    
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+M='Ge';
+if length(DB.textdata(1,:))==1  %% Octave data load
+
+    DB.textdata{1,1}=[DB.textdata{1,1} ',']; % patch, add a comma "," at the end
+    idxM=strfind(DB.textdata{1,1},',');
+    idx=strfind(DB.textdata{1,1},[',' M ',']);
+    idxM=find(idxM==idx);
+    
+    Ge = DB.data(:,idxM)';
+else  %% Matlab data load
+
+    for i=1:length(DB.textdata(1,:))
+      idx=strcmp(DB.textdata{1,i},M);
+      if idx==1
+        Ge=DB.data(:,i-1)';
+        % break % removing the break makes it slower but more compatible between Matlab and Octave
+      end
+    end 
+    
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+M='SiGe';
+if length(DB.textdata(1,:))==1  %% Octave data load
+
+    DB.textdata{1,1}=[DB.textdata{1,1} ',']; % patch, add a comma "," at the end
+    idxM=strfind(DB.textdata{1,1},',');
+    idx=strfind(DB.textdata{1,1},[',' M ',']);
+    idxM=find(idxM==idx);
+    
+    SiGe_bowing = DB.data(:,idxM)';
+else  %% Matlab data load
+
+    for i=1:length(DB.textdata(1,:))
+      idx=strcmp(DB.textdata{1,i},M);
+      if idx==1
+        SiGe_bowing=DB.data(:,i-1)';
+        % break % removing the break makes it slower but more compatible between Matlab and Octave
+      end
+    end 
+    
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+M='AlGaAs';
+if length(DB.textdata(1,:))==1  %% Octave data load
+
+    DB.textdata{1,1}=[DB.textdata{1,1} ',']; % patch, add a comma "," at the end
+    idxM=strfind(DB.textdata{1,1},',');
+    idx=strfind(DB.textdata{1,1},[',' M ',']);
+    idxM=find(idxM==idx);
+    
+    AlGaAs_bowing = DB.data(:,idxM)';
+else  %% Matlab data load
+
+    for i=1:length(DB.textdata(1,:))
+      idx=strcmp(DB.textdata{1,i},M);
+      if idx==1
+        AlGaAs_bowing=DB.data(:,i-1)';
+        % break % removing the break makes it slower but more compatible between Matlab and Octave
+      end
+    end 
+    
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+M='InGaAs';
+if length(DB.textdata(1,:))==1  %% Octave data load
+
+    DB.textdata{1,1}=[DB.textdata{1,1} ',']; % patch, add a comma "," at the end
+    idxM=strfind(DB.textdata{1,1},',');
+    idx=strfind(DB.textdata{1,1},[',' M ',']);
+    idxM=find(idxM==idx);
+    
+    InGaAs_bowing = DB.data(:,idxM)';
+else  %% Matlab data load
+
+    for i=1:length(DB.textdata(1,:))
+      idx=strcmp(DB.textdata{1,i},M);
+      if idx==1
+        InGaAs_bowing=DB.data(:,i-1)';
+        % break % removing the break makes it slower but more compatible between Matlab and Octave
+      end
+    end 
+    
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+M='AlInAs';
+if length(DB.textdata(1,:))==1  %% Octave data load
+
+    DB.textdata{1,1}=[DB.textdata{1,1} ',']; % patch, add a comma "," at the end
+    idxM=strfind(DB.textdata{1,1},',');
+    idx=strfind(DB.textdata{1,1},[',' M ',']);
+    idxM=find(idxM==idx);
+    
+    AlInAs_bowing = DB.data(:,idxM)';
+else  %% Matlab data load
+
+    for i=1:length(DB.textdata(1,:))
+      idx=strcmp(DB.textdata{1,i},M);
+      if idx==1
+        AlInAs_bowing=DB.data(:,i-1)';
+        % break % removing the break makes it slower but more compatible between Matlab and Octave
+      end
+    end 
+    
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+M='GaInP';
+if length(DB.textdata(1,:))==1  %% Octave data load
+
+    DB.textdata{1,1}=[DB.textdata{1,1} ',']; % patch, add a comma "," at the end
+    idxM=strfind(DB.textdata{1,1},',');
+    idx=strfind(DB.textdata{1,1},[',' M ',']);
+    idxM=find(idxM==idx);
+    
+    GaInP_bowing = DB.data(:,idxM)';
+else  %% Matlab data load
+
+    for i=1:length(DB.textdata(1,:))
+      idx=strcmp(DB.textdata{1,i},M);
+      if idx==1
+        GaInP_bowing=DB.data(:,i-1)';
+        % break % removing the break makes it slower but more compatible between Matlab and Octave
+      end
+    end 
+    
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+M='AlInP';
+if length(DB.textdata(1,:))==1  %% Octave data load
+
+    DB.textdata{1,1}=[DB.textdata{1,1} ',']; % patch, add a comma "," at the end
+    idxM=strfind(DB.textdata{1,1},',');
+    idx=strfind(DB.textdata{1,1},[',' M ',']);
+    idxM=find(idxM==idx);
+    
+    AlInP_bowing = DB.data(:,idxM)';
+else  %% Matlab data load
+
+    for i=1:length(DB.textdata(1,:))
+      idx=strcmp(DB.textdata{1,i},M);
+      if idx==1
+        AlInP_bowing=DB.data(:,i-1)';
+        % break % removing the break makes it slower but more compatible between Matlab and Octave
+      end
+    end 
+    
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+M='AlGaP';
+if length(DB.textdata(1,:))==1  %% Octave data load
+
+    DB.textdata{1,1}=[DB.textdata{1,1} ',']; % patch, add a comma "," at the end
+    idxM=strfind(DB.textdata{1,1},',');
+    idx=strfind(DB.textdata{1,1},[',' M ',']);
+    idxM=find(idxM==idx);
+    
+    AlGaP_bowing = DB.data(:,idxM)';
+else  %% Matlab data load
+
+    for i=1:length(DB.textdata(1,:))
+      idx=strcmp(DB.textdata{1,i},M);
+      if idx==1
+        AlGaP_bowing=DB.data(:,i-1)';
+        % break % removing the break makes it slower but more compatible between Matlab and Octave
+      end
+    end 
+    
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+M='GaInSb';
+if length(DB.textdata(1,:))==1  %% Octave data load
+
+    DB.textdata{1,1}=[DB.textdata{1,1} ',']; % patch, add a comma "," at the end
+    idxM=strfind(DB.textdata{1,1},',');
+    idx=strfind(DB.textdata{1,1},[',' M ',']);
+    idxM=find(idxM==idx);
+    
+    GaInSb_bowing = DB.data(:,idxM)';
+else  %% Matlab data load
+
+    for i=1:length(DB.textdata(1,:))
+      idx=strcmp(DB.textdata{1,i},M);
+      if idx==1
+        GaInSb_bowing=DB.data(:,i-1)';
+        % break % removing the break makes it slower but more compatible between Matlab and Octave
+      end
+    end 
+    
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+M='AlInSb';
+if length(DB.textdata(1,:))==1  %% Octave data load
+
+    DB.textdata{1,1}=[DB.textdata{1,1} ',']; % patch, add a comma "," at the end
+    idxM=strfind(DB.textdata{1,1},',');
+    idx=strfind(DB.textdata{1,1},[',' M ',']);
+    idxM=find(idxM==idx);
+    
+    AlInSb_bowing = DB.data(:,idxM)';
+else  %% Matlab data load
+
+    for i=1:length(DB.textdata(1,:))
+      idx=strcmp(DB.textdata{1,i},M);
+      if idx==1
+        AlInSb_bowing=DB.data(:,i-1)';
+        % break % removing the break makes it slower but more compatible between Matlab and Octave
+      end
+    end 
+    
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+M='AlGaSb';
+if length(DB.textdata(1,:))==1  %% Octave data load
+
+    DB.textdata{1,1}=[DB.textdata{1,1} ',']; % patch, add a comma "," at the end
+    idxM=strfind(DB.textdata{1,1},',');
+    idx=strfind(DB.textdata{1,1},[',' M ',']);
+    idxM=find(idxM==idx);
+    
+    AlGaSb_bowing = DB.data(:,idxM)';
+else  %% Matlab data load
+
+    for i=1:length(DB.textdata(1,:))
+      idx=strcmp(DB.textdata{1,i},M);
+      if idx==1
+        AlGaSb_bowing=DB.data(:,i-1)';
+        % break % removing the break makes it slower but more compatible between Matlab and Octave
+      end
+    end 
+    
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+M='GaAsSb';
+if length(DB.textdata(1,:))==1  %% Octave data load
+
+    DB.textdata{1,1}=[DB.textdata{1,1} ',']; % patch, add a comma "," at the end
+    idxM=strfind(DB.textdata{1,1},',');
+    idx=strfind(DB.textdata{1,1},[',' M ',']);
+    idxM=find(idxM==idx);
+    
+    GaAsSb_bowing = DB.data(:,idxM)';
+else  %% Matlab data load
+
+    for i=1:length(DB.textdata(1,:))
+      idx=strcmp(DB.textdata{1,i},M);
+      if idx==1
+        GaAsSb_bowing=DB.data(:,i-1)';
+        % break % removing the break makes it slower but more compatible between Matlab and Octave
+      end
+    end 
+    
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+M='InAsSb';
+if length(DB.textdata(1,:))==1  %% Octave data load
+
+    DB.textdata{1,1}=[DB.textdata{1,1} ',']; % patch, add a comma "," at the end
+    idxM=strfind(DB.textdata{1,1},',');
+    idx=strfind(DB.textdata{1,1},[',' M ',']);
+    idxM=find(idxM==idx);
+    
+    InAsSb_bowing = DB.data(:,idxM)';
+else  %% Matlab data load
+
+    for i=1:length(DB.textdata(1,:))
+      idx=strcmp(DB.textdata{1,i},M);
+      if idx==1
+        InAsSb_bowing=DB.data(:,i-1)';
+        % break % removing the break makes it slower but more compatible between Matlab and Octave
+      end
+    end 
+    
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+M='AlAsSb';
+if length(DB.textdata(1,:))==1  %% Octave data load
+
+    DB.textdata{1,1}=[DB.textdata{1,1} ',']; % patch, add a comma "," at the end
+    idxM=strfind(DB.textdata{1,1},',');
+    idx=strfind(DB.textdata{1,1},[',' M ',']);
+    idxM=find(idxM==idx);
+    
+    AlAsSb_bowing = DB.data(:,idxM)';
+else  %% Matlab data load
+
+    for i=1:length(DB.textdata(1,:))
+      idx=strcmp(DB.textdata{1,i},M);
+      if idx==1
+        AlAsSb_bowing=DB.data(:,i-1)';
+        % break % removing the break makes it slower but more compatible between Matlab and Octave
+      end
+    end 
+    
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+M='GaAsP';
+if length(DB.textdata(1,:))==1  %% Octave data load
+
+    DB.textdata{1,1}=[DB.textdata{1,1} ',']; % patch, add a comma "," at the end
+    idxM=strfind(DB.textdata{1,1},',');
+    idx=strfind(DB.textdata{1,1},[',' M ',']);
+    idxM=find(idxM==idx);
+    
+    GaAsP_bowing = DB.data(:,idxM)';
+else  %% Matlab data load
+
+    for i=1:length(DB.textdata(1,:))
+      idx=strcmp(DB.textdata{1,i},M);
+      if idx==1
+        GaAsP_bowing=DB.data(:,i-1)';
+        % break % removing the break makes it slower but more compatible between Matlab and Octave
+      end
+    end 
+    
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+M='InAsP';
+if length(DB.textdata(1,:))==1  %% Octave data load
+
+    DB.textdata{1,1}=[DB.textdata{1,1} ',']; % patch, add a comma "," at the end
+    idxM=strfind(DB.textdata{1,1},',');
+    idx=strfind(DB.textdata{1,1},[',' M ',']);
+    idxM=find(idxM==idx);
+    
+    InAsP_bowing = DB.data(:,idxM)';
+else  %% Matlab data load
+
+    for i=1:length(DB.textdata(1,:))
+      idx=strcmp(DB.textdata{1,i},M);
+      if idx==1
+        InAsP_bowing=DB.data(:,i-1)';
+        % break % removing the break makes it slower but more compatible between Matlab and Octave
+      end
+    end 
+    
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+M='AlAsP';
+if length(DB.textdata(1,:))==1  %% Octave data load
+
+    DB.textdata{1,1}=[DB.textdata{1,1} ',']; % patch, add a comma "," at the end
+    idxM=strfind(DB.textdata{1,1},',');
+    idx=strfind(DB.textdata{1,1},[',' M ',']);
+    idxM=find(idxM==idx);
+    
+    AlAsP_bowing = DB.data(:,idxM)';
+else  %% Matlab data load
+
+    for i=1:length(DB.textdata(1,:))
+      idx=strcmp(DB.textdata{1,i},M);
+      if idx==1
+        AlAsP_bowing=DB.data(:,i-1)';
+        % break % removing the break makes it slower but more compatible between Matlab and Octave
+      end
+    end 
+    
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+M='GaPSb';
+if length(DB.textdata(1,:))==1  %% Octave data load
+
+    DB.textdata{1,1}=[DB.textdata{1,1} ',']; % patch, add a comma "," at the end
+    idxM=strfind(DB.textdata{1,1},',');
+    idx=strfind(DB.textdata{1,1},[',' M ',']);
+    idxM=find(idxM==idx);
+    
+    GaPSb_bowing = DB.data(:,idxM)';
+else  %% Matlab data load
+
+    for i=1:length(DB.textdata(1,:))
+      idx=strcmp(DB.textdata{1,i},M);
+      if idx==1
+        GaPSb_bowing=DB.data(:,i-1)';
+        % break % removing the break makes it slower but more compatible between Matlab and Octave
+      end
+    end 
+    
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+M='InPSb';
+if length(DB.textdata(1,:))==1  %% Octave data load
+
+    DB.textdata{1,1}=[DB.textdata{1,1} ',']; % patch, add a comma "," at the end
+    idxM=strfind(DB.textdata{1,1},',');
+    idx=strfind(DB.textdata{1,1},[',' M ',']);
+    idxM=find(idxM==idx);
+    
+    InPSb_bowing = DB.data(:,idxM)';
+else  %% Matlab data load
+
+    for i=1:length(DB.textdata(1,:))
+      idx=strcmp(DB.textdata{1,i},M);
+      if idx==1
+        InPSb_bowing=DB.data(:,i-1)';
+        % break % removing the break makes it slower but more compatible between Matlab and Octave
+      end
+    end 
+    
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+M='AlPSb';
+if length(DB.textdata(1,:))==1  %% Octave data load
+
+    DB.textdata{1,1}=[DB.textdata{1,1} ',']; % patch, add a comma "," at the end
+    idxM=strfind(DB.textdata{1,1},',');
+    idx=strfind(DB.textdata{1,1},[',' M ',']);
+    idxM=find(idxM==idx);
+    
+    AlPSb_bowing = DB.data(:,idxM)';
+else  %% Matlab data load
+
+    for i=1:length(DB.textdata(1,:))
+      idx=strcmp(DB.textdata{1,i},M);
+      if idx==1
+        AlPSb_bowing=DB.data(:,i-1)';
+        % break % removing the break makes it slower but more compatible between Matlab and Octave
+      end
+    end 
+    
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\ No newline at end of file
diff --git a/PrintResults.m b/PrintResults.m
new file mode 100644
index 0000000..386c2d7
--- /dev/null
+++ b/PrintResults.m
@@ -0,0 +1,25 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Print Results %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%display('')
+display('===================================================')
+display('Intersubband Results:')
+display('===================================================')
+%display('')
+
+for i=1:length(Ec)
+  for j=1:length(Ec)
+    if j>i
+        display(strcat(...
+        'e',num2str(i),'-e',num2str(j),' = ',num2str( Ec(j)-Ec(i),'%.3f' ),'eV;   z'...
+        ,num2str(i),'-',num2str(j),' = ',num2str( z_dipole_c(i,j)*1e9,'%.3f' ),'nm;   f'...
+        ,num2str(i),'-',num2str(j),' = ',num2str( f_dipole_c(i,j),'%.3f' ) ...
+        )  )
+    end
+  end
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\ No newline at end of file
diff --git a/Results_GaAs-AlGaAs_0K.PNG b/Results_GaAs-AlGaAs_0K.PNG
new file mode 100644
index 0000000..f4fe3eb
Binary files /dev/null and b/Results_GaAs-AlGaAs_0K.PNG differ
diff --git a/Results_GaAs-AlGaAs_300K.PNG b/Results_GaAs-AlGaAs_300K.PNG
new file mode 100644
index 0000000..76e7645
Binary files /dev/null and b/Results_GaAs-AlGaAs_300K.PNG differ
diff --git a/Results_InGaAs-AlInAs_0K.PNG b/Results_InGaAs-AlInAs_0K.PNG
new file mode 100644
index 0000000..0390f71
Binary files /dev/null and b/Results_InGaAs-AlInAs_0K.PNG differ
diff --git a/Results_InGaAs-AlInAs_300K.PNG b/Results_InGaAs-AlInAs_300K.PNG
new file mode 100644
index 0000000..62b543a
Binary files /dev/null and b/Results_InGaAs-AlInAs_300K.PNG differ
diff --git a/Results_InGaAs-GaAs_300K.PNG b/Results_InGaAs-GaAs_300K.PNG
new file mode 100644
index 0000000..580bbc0
Binary files /dev/null and b/Results_InGaAs-GaAs_300K.PNG differ
diff --git a/Schrod_2bands_Kane_f.m b/Schrod_2bands_Kane_f.m
new file mode 100644
index 0000000..779fb25
--- /dev/null
+++ b/Schrod_2bands_Kane_f.m
@@ -0,0 +1,121 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%% Schrodinger solver on uniform grid with m(z,E)!!! %%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%% With the non-parabolic band 2x2k.p Kane model %%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function[Ec,psi_c]=Schrod_2bands_Kane_f(z,Vc,Eg,EP,Dso,n,ac,av,bv,exx,ezz)
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Constants %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+h=6.62606896E-34;               %% Planck constant J.s
+hbar=h/(2*pi);
+e=1.602176487E-19;              %% charge de l electron Coulomb
+m0=9.10938188E-31;              %% electron mass kg
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+Nz=length(z);
+dz = z(2)-z(1);
+
+eyy = exx;
+DCBO   = -abs(ac).*(exx+eyy+ezz) ; % shift of the CB due to strain
+DVBOLH = +abs(av).*(exx+eyy+ezz) + abs(bv).*(exx-ezz) ; % shift of the VB due to strain
+DVBOSO = +abs(av).*(exx+eyy+ezz) ; % shift of the VB due to strain
+
+Vc=Vc+DCBO;
+Vc(1)=5;
+Vc(end)=5;
+
+shift=min(Vc);
+Vc=Vc-shift;
+
+Vv=Vc-Eg+DVBOLH;
+Vso=Vc-Dso-Eg+DVBOSO;
+Vveff= (2*Vv+Vso)/3;  % here is an effective valence band that make a ratio between LH and SO
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%% Building of the operators %%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%% First derivative %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%DZ1c = (0.5)*diag(ones(1,Nz-1),+1) + (-0.5)*diag(ones(1,Nz-1),-1)  ;
+DZ1b =   (1)*diag(ones(1,Nz)  ,0 ) +   (-1)*diag(ones(1,Nz-1),-1)  ;
+DZ1f =   (1)*diag(ones(1,Nz-1),+1) +   (-1)*diag(ones(1,Nz)  ,0 )  ;
+
+%DZ1c=DZ1c/dz;
+DZ1b=DZ1b/dz;
+DZ1f=DZ1f/dz;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%% Second derivative %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+DZ2 =(-2)*diag(ones(1,Nz)) + (1)*diag(ones(1,Nz-1),-1) + (1)*diag(ones(1,Nz-1),1);
+DZ2=DZ2/dz^2;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%% Building of the Hamiltonien %%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+Vc    = [ (Vc(1:end-1)    + Vc(2:end))   / 2   Vc(end)    ];
+Vveff = [ (Vveff(1:end-1) + Vveff(2:end))/ 2   Vveff(end) ];
+EP    = [ (EP(1:end-1)    + EP(2:end))   / 2   EP(end)    ];
+
+H0=(-(hbar^2)/(2*m0)) *  DZ2   ;
+
+H11 =  H0 + diag(Vc*e) ;
+H22 = -H0 + diag( (Vveff)*e ) ;
+
+% Xunpeng Ma et al. JAP, 114, 063101 (2013)
+% "Two-band finite difference method for the bandstructure calculation with nonparabolicity effects in quantum cascade lasers"
+% ==> It s seems to be by far the most accurate method
+
+H12 = +1*hbar/sqrt(2*m0) * (  diag(sqrt(EP*e),0) + diag(sqrt(EP(1:end-1)*e),-1)  ) .* DZ1b ;
+H21 = -1*hbar/sqrt(2*m0) * (  diag(sqrt(EP*e),0) + diag(sqrt(EP(1:end-1)*e),+1)  ) .* DZ1f ;
+
+H2x2=[
+H11   H12
+H21   H22
+];
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%% Diagonalisation of the Hamiltonien %%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+H2x2=sparse(H2x2);
+[psi_2x2,Energy] = eigs(H2x2,n,'SM');
+E_2x2 = diag(Energy)/e ;
+
+psi_c=[];
+Ec=[];
+
+for i=1:n
+  if E_2x2(i) > min(Vc)
+    Ec=[Ec abs(E_2x2(i))];
+    psi_c = [psi_c  psi_2x2(1:length(z),i) ];
+  end
+end
+
+Ec=Ec'+shift;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%% Normalization of the Wavefunction %%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+for i=1:length(Ec)
+    psi_c(:,i)=psi_c(:,i)/sqrt(trapz(z',abs(psi_c(:,i)).^2)); 
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% here is a small patch due to differences between Octave and Matlab
+% Matlab order the eigen values while Octave reverse it
+
+if length(Ec)>1
+if Ec(1)>Ec(2)
+  psi_c=psi_c(:,end:-1:1);
+  Ec=Ec(end:-1:1);
+end
+end
+
+end
\ No newline at end of file
diff --git a/Schrod_Nbands_shoot_f.m b/Schrod_Nbands_shoot_f.m
new file mode 100644
index 0000000..55988dd
--- /dev/null
+++ b/Schrod_Nbands_shoot_f.m
@@ -0,0 +1,116 @@
+function[E,psi]=Schrod_Nbands_shoot_f(z,V0,me,n,Evec,dE,precision)
+
+method=2;         % method "2" is much more acurate than "1" but take also more time...
+e=min(V0);
+Emax=max(V0)+0.1;
+
+C=0; N=0; E=[];
+psi=z*0+1; psi_old=psi;
+
+
+while e<Emax && length(E)<n
+    
+    C = C+1;
+    
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% In this block, it is scanning to find the mass at the right Energy, m(E)
+
+    idx=[0 0];
+    epsi=1;
+    while length(idx)>1
+        idx=find( abs( (Evec-e )) < epsi ) ;
+      if isempty(idx)
+          idx=IDX(1);
+        else
+          IDX=idx;
+          epsi=epsi/2;
+        end
+    end
+    Mass=me(idx,:);
+%    [ e Mass(round(end/2)) ];
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    
+    psi_old = psi;
+    psi=Schroed1D_Euler_Eval(z,V0,Mass,e,method);
+    
+    if (sign( psi(end) ) ~= sign( psi_old(end) ) ) && C>1 % here, I catch a quantum state because the last point of psi change sign
+    
+        N=N+1; de=dE;
+        
+        while abs(de)>precision
+            if sign( psi(end) ) ~= sign( psi_old(end) )
+                 de = -de/2;
+            end
+            e=e+de;
+            psi_old=psi;
+            psi=Schroed1D_Euler_Eval(z,V0,Mass,e,method);
+        end
+
+        E(N,:)=e; PSI(:,N)= psi;
+        C=0;
+    end	
+    
+    e=e+dE;
+end
+
+psi=PSI;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%% Normalization of the Wavefunction %%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+for i=1:N
+    z=z';
+    psi(:,i)=psi(:,i)/sqrt(trapz(z,abs(psi(:,i)).^2));  
+end
+
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+function[psi]=Schroed1D_Euler_Eval(z,V0,Mass,ee,method)
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+h=6.62606896E-34;               %% Planck constant [J.s]
+hbar=h/(2*pi);
+e=1.602176487E-19;              %% electron charge [C]
+m0=9.10938188E-31;              %% electron mass [kg]
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+dz=z(2)-z(1);
+dz=[ dz diff(z)];
+y = [0 1];
+
+if method==1
+    
+    for i = 1:length(z)                                  %% number of z steps to take
+        dy(2) = -2*e/(hbar^2) * ( ee  - V0(i) ) * y(1) ; %% Equation for dv/dz
+        dy(1) =  Mass(i)*m0*y(2);                        %% Equation for dx/dz
+        y = y + dz(i)*dy ;                               %% integrate both equations with Euler
+        psi(i)= y(1);
+    end
+    
+elseif method==2
+    
+    for i = 1:length(z)                                  %% number of z steps to take
+    
+        dy(2) = -2*e/(hbar^2) * ( ee  - V0(i) ) * y(1) ; %% Equation for dv/dz
+        dy(1) =  Mass(i)*m0*y(2);                        %% Equation for dx/dz
+        
+        K = y + 0.5*dz(i)*dy;
+        dK(2) = -2*e/(hbar^2) * ( ee  - V0(i) ) * K(1) ; %% Equation for dv/dz
+        dK(1) =  Mass(i)*m0*K(2);                        %% Equation for dx/dz
+        
+        y =  y + dz(i)*dK;
+        psi(i)= y(1);
+    end
+end
+
+end
\ No newline at end of file
diff --git a/SchrodingerPoisson1D_CB_Kane_Main.m b/SchrodingerPoisson1D_CB_Kane_Main.m
new file mode 100644
index 0000000..68b3466
--- /dev/null
+++ b/SchrodingerPoisson1D_CB_Kane_Main.m
@@ -0,0 +1,491 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%% last update 20December2021, lne %%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% This program solves the Schrodinger-Poisson equation in the conduction band 
+% for any heterostructures in the Zinc-Blend material.
+% The program solves the Schrodinger equation with m(E,z) using different algorithms:
+% -> FDM (Finite Difference Method) 
+%          and / or
+% -> Scaning/Shooting method
+% The non-parabolicity is implemented via the kp model using the Kane model
+% A strain model is included. It basically shifts the conduction band edge
+% The strain is mainly interesting for InGaAs/GaAs heterostructures
+% The non-parabolicity is also included into the density of states for the Poisson solver.
+% It follows the book of Paul Harrison and actually makes the 2d density of states not constant
+% -> Additionnal material can be added in the "materialDB_ZB.csv" file
+% -> II-VI and cubic nitride material parameters are available but should
+% be grabt in the "Library.m" file
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% If the code doesn t converge:
+% -> decrease the doping
+% -> increase the resolution dz
+% -> increase the temperature (T=0K is very bad while T=10K is already much better)
+% -> increase the amount of loops, Nloops
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+clear all
+close all
+clc
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Constants %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+h    = 6.62606896E-34;              %% Planck constant [J.s]
+hbar = h/(2*pi);
+e    = 1.602176487E-19;             %% electron charge [C]
+m0   = 9.10938188E-31;              %% electron mass [kg]
+Epsi0= 8.854187817620E-12;          %% Vaccum dielectric constant [F/m]
+kB   = 1.3806488E-23;               %% Boltzmann's constant [J/K]
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+Nloops = 5;                   % number of loops
+model  = 2;                   % 1=shooting; 2=Kane
+StrainModel = 0;              % Activate Strain model
+n      = 3;                   % number of solution asked per model
+ScF    = 0.1;                 % scaling factor to plot the wave function [Without Dimension]
+dz     = 1e-10;               % resolution of the grid [m]
+T      = 300;                 % Temperature [Kelvin]
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+DisplayResults   = 1;         % Switch to print or not the ISB dipoles on the shell
+
+plot_density     = 1;         % Activate the plot 0 or 1
+plot_convergence = 0;         % Activate the plot 0 or 1
+plot_field       = 0;         % Activate the plot 0 or 1
+plot_Vbending    = 0;         % Activate the plot 0 or 1
+plot_mass        = 0;         % Activate the plot 0 or 1
+plot_ro          = 0;         % Activate the plot 0 or 1
+plot_Epsi        = 0;         % Activate the plot 0 or 1
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+Library;                  % load material parameter DB from "materialDB_ZB.csv"
+ExtractParameters;        % extract parameter from the Library
+TernaryAlloy;             % compute the ternary alloy    
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%% import the layer structure file %%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% first column is the material used from the "library"
+% second column is the length of the layer in nm
+% third column is the n doping volumique of that layer in 1e18 cm-3
+
+% You have to put a resonable amount of doping! Otherwise, it will diverge 
+% and you will have to damp it more by increasing the number of loops !
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+input_file;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%% NOTHING TO CHANGE ANYMORE !!! %%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%% Grabbing the parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+zt   = M(:,end-1)*1e-9;     % conversion of the length from nm to meter
+Dopt = M(:,end)*1e18*1e6;   % n doping conversion from cm-3 to m-3
+
+Egt  = M(:,idx_Eg6c) - (M(:,idx_alphaG)*T^2) ./ (T+M(:,idx_betaG));   %Eg = Eg0 - (a*T.^2)./(T + b);
+CBOt = Egt+M(:,idx_VBO);	  % CBO form band gap difference and temperature
+Dsot = M(:,idx_Dso);        % Spin-Orbit shift band parameter
+EPt_K= M(:,idx_EP_K);       % EP Kane
+Epsit= M(:,idx_Epsi);       % Epsilon; dielectric constant
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Strain Model %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+at   = M(:,idx_a);           % lattice parameter
+act  = M(:,idx_ac);          % Conduction band strain offset parameter
+avt  = M(:,idx_av);          % Valence band strain offset parameter
+bvt  = M(:,idx_bv);          % Valence band strain offset parameter
+c11t = M(:,idx_c11);        % strain parameter
+c12t = M(:,idx_c12);        % strain parameter
+
+a0   = substrate(idx_a);
+
+if StrainModel == 1
+  exxt =  (a0-at)/a0; % eyyt =  exxt;
+  ezzt = -2*c12t./c11t.*exxt;
+else
+  exxt =  (a0-at)/a0 * 0; % eyyt =  exxt;
+  ezzt = -2*c12t./c11t.*exxt;
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Discretisation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% here, I descretize the grid z, the potential V0 and that values that are needed
+
+z=0; V0=CBOt(1); Eg=Egt(1); Dso=Dsot(1); EP_K=EPt_K(1);
+Dop=Dopt(1); Epsi=Epsit(1);
+ac=act(1); av=avt(1); bv=bvt(1); exx=exxt(1); ezz=ezzt(1);
+
+for i=1:length(zt)
+    t=zt(i);
+    zv= (z(end)+dz): dz : (z(end)+dz)+t;
+    z=[z zv];
+    V0  = [ V0     ones(size(zv)) * CBOt(i)  ];
+    Eg  = [ Eg     ones(size(zv)) * Egt(i)   ];
+    Dop = [ Dop    ones(size(zv)) * Dopt(i)  ];
+    Epsi= [ Epsi   ones(size(zv)) * Epsit(i) ];
+    EP_K= [ EP_K   ones(size(zv)) * EPt_K(i) ];
+    Dso = [ Dso    ones(size(zv)) * Dsot(i)  ];
+    ac  = [ ac     ones(size(zv)) * act(i)   ];
+    av  = [ av     ones(size(zv)) * avt(i)   ];
+    bv  = [ bv     ones(size(zv)) * bvt(i)   ];
+    exx = [ exx    ones(size(zv)) * exxt(i)  ];
+    ezz = [ ezz    ones(size(zv)) * ezzt(i)  ];
+end
+
+V0=V0-min(V0);             % Shift the band in order to get the bottom of the well at zero
+V0=(F0*z)+V0;              % adding the electric field to the potential
+
+Ltot=z(end)-z(1);
+
+eyy  = exx;
+DCBO = -abs(ac).*(exx+eyy+ezz) ;                      % shift of the CB due to strain
+
+Ntott= Dopt.*zt;
+Ntot = sum(Ntott);   % total number of charges
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+dE1=1e-3;
+dE2=1e-2;
+    
+if StrainModel == 0
+    E1 = min(V0):dE1:min(V0)+0.2 ;
+    E2 = E1(end):dE2:max(V0);
+    En = sort([E1 E2]);
+else
+    E1 = min(V0+DCBO):dE1:min(V0+DCBO)+0.2 ;
+    E2 = E1(end):dE2:max(V0+DCBO);
+    En = sort([E1 E2]);
+end
+
+if model==1
+    EEn     = repmat(En',  [1 length(z)]);
+    V0mat   = repmat(V0 ,  [length(En) 1]);
+    Egmat   = repmat(Eg ,  [length(En) 1]);
+    Dsomat  = repmat(Dso,  [length(En) 1]);
+    EPmat   = repmat(EP_K ,[length(En) 1]);
+    meK2    = 1 ./ ( 1 +  EPmat ./ (EEn + Egmat-V0mat+1/3*Dsomat)  ); % m(z,E) in the 2 bands Kane s model
+    
+    % here, I linearize the mass versus the energy
+    meff0   = 1 ./ ( 1 +  EPmat ./ (Egmat-V0mat+1/3*Dsomat)  ); 
+    alpha   = ( meK2(end,:)./meff0(end,:) - 1 ) ./ EEn(end,:);
+    alphamat= repmat(alpha , [length(En) 1] );
+    meK2lin = meff0.*(1+alphamat.*EEn); % => linearized mass
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%% Starting of the Poisson s loop %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+Vs=zeros(size(z)); Vsold=Vs;
+ntot=0; nloop=1;
+ErrVec=1; sumVtotVec=1;
+
+while nloop<Nloops
+    
+    nloop
+    x = 0.5;
+    Vbending=Vs*x + Vsold*(1-x);
+    Vtot=V0+Vbending;
+        
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%% schrodinger solver %%%%%%%%%%%%%%%%%%%%%%%%%%%%
+        
+    if model==1
+        dE=0.002;
+        precision=1e-7;
+        Vtotmat = repmat(Vtot, [length(En) 1]);
+        meK2    = 1 ./ ( 1 +  EPmat ./ (EEn + Egmat-Vtotmat+1/3*Dsomat)  ); % m(z,E) in the 2 bands Kane s model
+        [Ec,psic] = Schrod_Nbands_shoot_f(z,Vtot+DCBO,meK2,n,En,dE,precision);
+    elseif model==2
+        [Ec,psic] = Schrod_2bands_Kane_f(z,Vtot,Eg,EP_K,Dso,n,ac,av,bv,exx,ezz);
+    end
+
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    % Here, I re-define the energy grid in order optimize the meshing
+    dE1=1e-4;
+    dE2=1e-2;
+    if StrainModel == 0
+        E1 = Ec(1):dE1:Ec(1)+0.1 ;
+        E2 = E1(end):dE2:max(Vtot);
+        En=sort([E1 E2]);
+    else
+        E1 = Ec(1):dE1:Ec(1)+0.1 ;
+        E2 = E1(end):dE2:max(Vtot+DCBO);
+        En=sort([E1 E2]);
+    end
+
+    EEn     = repmat(En' , [1 length(z) ]);
+    V0mat   = repmat(V0  , [length(En) 1]);
+    Vtotmat = repmat(Vtot, [length(En) 1]);
+    Egmat   = repmat(Eg  , [length(En) 1]);
+    Dsomat  = repmat(Dso , [length(En) 1]);
+    EPmat   = repmat(EP_K, [length(En) 1]);
+    meK2    = 1 ./ ( 1 +  EPmat ./ (EEn + Egmat-Vtotmat+1/3*Dsomat)  ); % m(z,E) in the 2 bands Kane s model
+    
+    % here, I linearize the mass versus the energy
+    meff0   = 1 ./ ( 1 +  EPmat ./ (Egmat-Vtotmat+1/3*Dsomat)  ); 
+    alpha   = ( meK2(end,:)./meff0(end,:) - 1 ) ./ EEn(end,:);
+    alphamat= repmat(alpha , [length(En) 1] );
+    meK2lin = meff0.*(1+alphamat.*EEn); % => linearized mass
+    
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    % Paul Harrisson
+    % Quantum Wells, Wires and Dots.
+    % 4th edition (2016),
+    % chap 2 : "Solutions to Schrodinger's equation"
+    % 2.42: "Two-dimensional systems" page 31
+    
+    meffro = meff0.*(1+2*alphamat.*EEn); % => linearized mass
+    ro=[];
+    for i=1:length(Ec)
+        %ro( En>Ec(i),:,i) = e*   meK2(En>Ec(i),:) * m0/(pi*hbar^2);
+        %ro( En>Ec(i),:,i) = e*meK2lin(En>Ec(i),:) * m0/(pi*hbar^2);
+        ro( En>Ec(i),:,i) = e*meffro(En>Ec(i),:) * m0/(pi*hbar^2);
+        ro( En<Ec(i),:,i) = 0;
+    end
+    
+    % here is what make the code very slow
+    if Ntot==0
+        idx=find(min(CBOt)==CBOt);
+        Ef=-Egt(idx(1))/2;
+        NN=0*Ec';
+        roEf=ro*0;
+    else
+        [Ef,NN,roEf]=find_Ef_f(z,Ec,psic,En,ro,Ntot,T);
+    end
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    
+    ntot2 = repmat(NN,[length(z) 1]).*abs(psic).^2 ;
+    
+    ntot = sum(ntot2,2)' - Dop;  % remove the charge positives (ions)
+        
+    %%%%%%%%%%%%%%%%%%%%%%%%%% Electrical Field %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    F = e*cumtrapz(z,ntot)./(Epsi0*Epsi);
+    MF  = trapz(z,F)/(z(end)-z(1));  % MF=mean(F) on a nonlinear grid z
+    F = F-MF;
+    %%%%%%%%%%%%%%%%%%%%%%%%%%% New Potential %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    Vsold = Vs;
+    Vs    = -cumsum(F)*dz;
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%% Convergence %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+    Err = abs(  1 - sumVtotVec(end)/sum(Vs)  );
+    sumVtotVec(nloop) = sum(Vs);
+    ErrVec = [ErrVec Err];
+
+    nloop=nloop+1;
+
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Display Results %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+scale_PSI;
+computesISBdipoles;
+
+if DisplayResults == 1
+  PrintResults;
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figures %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%figure('position',[-3500 100 1000 700],'color','w');
+figure('position',[10 100 1000 700],'color','w');
+subplot(1,1,1,'fontsize',15)
+hold on; grid on; box on;
+col=colormap(jet);
+
+if plot_density==1
+    grid off
+    pcolor(z*1e9,En,sum(ROEf,3)*1e-6 )
+    set(gca,'color',col(1,:))
+    shading flat
+    hcb=colorbar;
+    caxis([0 max(max(sum(ROEf,3)*1e-6))+1])
+    title(hcb,'\fontsize{8}cm-3')
+    
+    if StrainModel == 0
+        plot(z*1e9,V0,  'w--','linewidth',1)
+        plot(z*1e9,Vtot,'w-' ,'linewidth',1)
+    else
+        plot(z*1e9,V0+DCBO,  'k--','linewidth',1)
+        plot(z*1e9,Vtot+DCBO,'k-' ,'linewidth',1)
+    end
+elseif plot_density==0
+    if StrainModel == 0
+        plot(z*1e9,V0,  'b--','linewidth',1)
+        plot(z*1e9,Vtot,'b-' ,'linewidth',1)
+    else
+        plot(z*1e9,V0+DCBO,  'k--','linewidth',1)
+        plot(z*1e9,Vtot+DCBO,'k-' ,'linewidth',1)
+    end
+end
+
+for i=1:length(Ec)
+    plot(z*1e9,PSIc(:,i),'color','r','linewidth',1)
+end
+
+plot([z(1) z(end)]*1e9,[1 1]*Ef,'g','linewidth',1)
+text(z(end)*1e9*0.95,Ef+0.01,'\color{green}Fermi')
+
+xlabel('z (nm)')
+ylabel('Energy (eV)')
+xlim([0 z(end)*1e9])
+if StrainModel == 1
+  title(strcat('T=',num2str(T),'K ; dz=',num2str(dz*1e9),'nm; Ntot=',num2str(Ntot*1e-4,'%.1e'),'cm-2; with STRAIN'))
+else
+  title(strcat('T=',num2str(T),'K ; dz=',num2str(dz*1e9),'nm; Ntot=',num2str(Ntot*1e-4,'%.1e'),'cm-2; without STRAIN'))
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+if plot_convergence==1
+    figure('color','w')
+    semilogy(1:nloop,ErrVec,'bo-')
+    hold on; grid on; box on;
+    xlabel('Cycles')
+    ylabel('Convergence (norm. units)')
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+if plot_field==1
+    
+    figure('color','w')
+    hold on; grid on; box on;
+    [AX,H1,H2]=plotyy(z*1e9,F*1e-2*1e-3,z*1e9,Dop*1e-18*1e-6);
+        
+    set(H1,'color','r')
+    set(H2,'color','b')
+    
+    xlabel('z (nm)')
+    ylabel(AX(1),'E- field (kV/cm)','color','red')
+    ylabel(AX(2),'Doping (1e18 cm-3)','color','blue')
+    
+    set(AX(1),'ycolor','red')
+    set(AX(2),'ycolor','blue')
+        
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+if plot_Vbending==1
+
+    figure('color','w')
+    hold on; grid on; box on;
+
+    [AX,H1,H2]=plotyy(z*1e9,Vbending,z*1e9,ntot*1e-18*1e-6);
+
+    set(H1,'color','r')
+    set(H2,'color','b')
+    
+    xlabel('z (nm)')
+    ylabel(AX(1),'Vbending (eV)','color','red')
+    ylabel(AX(2),'ntot (1e18 cm-3)','color','blue')
+    
+    set(AX(1),'ycolor','red')
+    set(AX(2),'ycolor','blue')
+
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+if plot_mass==1
+    figure('color','w')
+    
+    subplot(1,2,1)
+    hold on; grid on; box on;
+    pcolor(z*1e9,En,meK2)
+    
+    colormap(jet)
+    hcb=colorbar;
+    title(hcb,'\fontsize{8}meff')
+    shading flat
+    xlabel('z (nm)')
+    ylabel('Energy (eV)')
+    title('meK2')
+    xlim([0 z(end)*1e9])
+    
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    
+    subplot(1,2,2)
+    hold on; grid on; box on;
+    pcolor(z*1e9,En,meK2lin)
+    
+    colormap(jet)
+    hcb=colorbar;
+    title(hcb,'\fontsize{8}meff')
+    shading flat
+    xlabel('z (nm)')
+    ylabel('Energy (eV)')
+    title('meK2-lin')
+    xlim([0 z(end)*1e9])
+    
+%     figure
+%     hold on;grid on;
+%     idx=250;
+%     plot(En,meK2(:,idx),'b')
+%     plot(En,meK2lin(:,idx),'r')
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+if plot_ro==1
+    figure('color','w')
+    hold on; grid on; box on;
+    surf(z*1e9,En,sum(ro,3))
+    view(45,30)
+    colormap(jet)
+    hcb=colorbar;
+    title(hcb,'\fontsize{8}ro')
+    shading flat
+    xlabel('z (nm)')
+    ylabel('Energy (eV)')
+    xlim([0 z(end)*1e9])
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+if plot_Epsi==1
+    figure('color','w')
+    hold on; grid on; box on;
+    plot(z*1e9,Epsi,'.-')
+    xlabel('z (nm)')
+    ylabel('Epsilon')
+    xlim([0 z(end)*1e9])
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/SchrodingerPoisson1D_CB_Main.m b/SchrodingerPoisson1D_CB_Main.m
new file mode 100644
index 0000000..19745e8
--- /dev/null
+++ b/SchrodingerPoisson1D_CB_Main.m
@@ -0,0 +1,456 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%% last update 25December2021, lne %%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% This program solves the Schrodinger-Poisson equation in the conduction band 
+% for any heterostructures.
+% The program solves the Schrodinger equation with m(E,z) using the Scaning/Shooting
+% method. The non-parabolicity is implemented via the alpha parameter, alpha=1/Egap;
+% meff(E)=meff(0)*(1+alpha*E); 
+% It can be easily modified by tuning the alpha in each layer => line 109
+% A strain model is included. It basically shifts the conduction band edge
+% The strain is mainly interesting for InGaAs/GaAs heterostructures
+% The non-parabolicity is also included into the density of states for the Poisson solver.
+% It follows the book of Paul Harrison and actually makes the 2d density of states not constant
+% -> Additionnal material can be added in the "materialDB_ZB.csv" file
+% -> II-VI and cubic nitride material parameters are available but should
+% be grabt in the "Library.m" file
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% If the code doesn t converge:
+% -> decrease the doping
+% -> increase the resolution dz
+% -> increase the temperature (T=0K is very bad while T=10K is already much better)
+% -> increase the amount of loops, Nloops
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+clear all
+close all
+clc
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Constants %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+h    = 6.62606896E-34;              %% Planck constant [J.s]
+hbar = h/(2*pi);
+e    = 1.602176487E-19;             %% electron charge [C]
+m0   = 9.10938188E-31;              %% electron mass [kg]
+Epsi0= 8.854187817620E-12;          %% Vaccum dielectric constant [F/m]
+kB   = 1.3806488E-23;               %% Boltzmann's constant [J/K]
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+Nloops = 5;                   % number of loops
+StrainModel = 0;              % Activate Strain model
+n      = 4;                   % number of solution asked per model
+ScF    = 0.1;                 % scaling factor to plot the wave function [Without Dimension]
+dz     = 1e-10;               % resolution of the grid [m]
+T      = 300;                 % Temperature [Kelvin]
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+DisplayResults   = 1;         % Switch to print or not the ISB dipoles on the shell
+
+plot_density     = 1;         % Activate the plot 0 or 1
+plot_convergence = 0;         % Activate the plot 0 or 1
+plot_field       = 0;         % Activate the plot 0 or 1
+plot_Vbending    = 0;         % Activate the plot 0 or 1
+plot_mass        = 0;         % Activate the plot 0 or 1
+plot_ro          = 0;         % Activate the plot 0 or 1
+plot_Epsi        = 0;         % Activate the plot 0 or 1
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+Library;                  % load material parameter DB from "materialDB_ZB.csv"
+ExtractParameters;        % extract parameter from the Library
+TernaryAlloy;             % compute the ternary alloy    
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%% import the layer structure file %%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% first column is the material used from the "library"
+% second column is the length of the layer in nm
+% third column is the n doping volumique of that layer in 1e18 cm-3
+
+% You have to put a resonable amount of doping! Otherwise, it will diverge 
+% and you will have to damp it more by increasing the number of loops !
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+input_file;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%% NOTHING TO CHANGE ANYMORE !!! %%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%% Grabbing the parameters %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+zt    = M(:,end-1)*1e-9;      % conversion of the length from nm to meter
+Dopt  = M(:,end)*1e18*1e6;    % n doping conversion from cm-3 to m-3
+
+Egt   = M(:,idx_Eg6c) - (M(:,idx_alphaG)*T^2) ./ (T+M(:,idx_betaG));   %Eg = Eg0 - (a*T.^2)./(T + b);
+CBOt  = Egt+M(:,idx_VBO);     % CBO form band gap difference and temperature
+%Dsot = M(:,idx_Dso);        % Spin-Orbit shift band parameter
+%EPt_K= M(:,idx_EP_K);       % EP Kane
+Epsit = M(:,idx_Epsi);        % Epsilon; dielectric constant
+
+mefft = M(:,idx_me);          % electron effective mass at band edge
+alphat= 1./Egt;               % non-parabolicity parameter => can be adjusted!
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Strain Model %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+at  = M(:,idx_a);           % lattice parameter
+act = M(:,idx_ac);          % Conduction band strain offset parameter
+avt = M(:,idx_av);          % Valence band strain offset parameter
+bvt = M(:,idx_bv);          % Valence band strain offset parameter
+c11t = M(:,idx_c11);        % strain parameter
+c12t = M(:,idx_c12);        % strain parameter
+
+a0   = substrate(idx_a);
+
+if StrainModel == 1
+  exxt =  (a0-at)/a0; % eyyt =  exxt;
+  ezzt = -2*c12t./c11t.*exxt;
+else
+  exxt =  (a0-at)/a0 * 0; % eyyt =  exxt;
+  ezzt = -2*c12t./c11t.*exxt;
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Discretisation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% here, I descretize the grid z, the potential V0 and that values that are needed
+
+z=0; V0=CBOt(1); Eg=Egt(1); 
+meff=mefft(1); alpha=alphat(1);
+Dop=Dopt(1); Epsi=Epsit(1);
+ac=act(1); av=avt(1); bv=bvt(1); exx=exxt(1); ezz=ezzt(1);
+
+for i=1:length(zt)
+    t    = zt(i);
+    zv   = (z(end)+dz): dz : (z(end)+dz)+t;
+    z    = [z zv];
+    V0   = [ V0     ones(size(zv)) * CBOt(i)  ];
+    Eg   = [ Eg     ones(size(zv)) * Egt(i)   ];
+    Dop  = [ Dop    ones(size(zv)) * Dopt(i)  ];
+    Epsi = [ Epsi   ones(size(zv)) * Epsit(i) ];
+    meff = [ meff   ones(size(zv)) * mefft(i) ];
+    alpha= [ alpha  ones(size(zv)) * alphat(i)];
+    ac   = [ ac     ones(size(zv)) * act(i)   ];
+    av   = [ av     ones(size(zv)) * avt(i)   ];
+    bv   = [ bv     ones(size(zv)) * bvt(i)   ];
+    exx  = [ exx    ones(size(zv)) * exxt(i)  ];
+    ezz  = [ ezz    ones(size(zv)) * ezzt(i)  ];
+end
+
+V0=V0-min(V0);             % Shift the band in order to get the bottom of the well at zero
+V0=(F0*z)+V0;              % adding the electric field to the potential
+
+Ltot=z(end)-z(1);
+
+eyy = exx;
+DCBO   = -abs(ac).*(exx+eyy+ezz) ;                      % shift of the CB due to strain
+
+Ntott=Dopt.*zt;
+Ntot=sum(Ntott);   % total number of charges
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+dE1=1e-3;
+dE2=1e-2;
+    
+if StrainModel == 0
+    E1 = min(V0):dE1:min(V0)+0.2 ;
+    E2 = E1(end):dE2:max(V0);
+    En=sort([E1 E2]);
+else
+    E1 = min(V0+DCBO):dE1:min(V0+DCBO)+0.2 ;
+    E2 = E1(end):dE2:max(V0+DCBO);
+    En=sort([E1 E2]);
+end
+
+EEn     = repmat(En',   [1  length(z)]);
+V0mat   = repmat(V0 ,   [length(En) 1]);
+Egmat   = repmat(Eg ,   [length(En) 1]);
+meffmat = repmat(meff , [length(En) 1]);
+alphamat= repmat(alpha, [length(En) 1]);
+melin   = meffmat.*(1+alphamat.*(EEn-V0mat)); % => linearized mass
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%% Starting of the Poisson s loop %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+Vs=zeros(size(z)); Vsold=Vs;
+ntot=0; nloop=1;
+ErrVec=1; sumVtotVec=1;
+
+while nloop<Nloops
+    
+    nloop
+    x = 0.5;
+    Vbending=Vs*x + Vsold*(1-x);
+    Vtot=V0+Vbending;
+        
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%% schrodinger solver %%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    
+    dE=0.002; precision=1e-7;
+    Vtotmat = repmat(Vtot, [length(En) 1]);
+    melin = meffmat.*(1+alphamat.*(EEn-Vtotmat)); % => linearized mass
+    [Ec,psic] = Schrod_Nbands_shoot_f(z,Vtot+DCBO,melin,n,En,dE,precision);
+
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    % Here, I re-define the energy grid in order optimize the meshing
+    dE1=1e-4;
+    dE2=1e-2;
+    if StrainModel == 0
+        E1 = Ec(1):dE1:Ec(1)+0.1 ;
+        E2 = E1(end):dE2:max(Vtot);
+        En=sort([E1 E2]);
+    else
+        E1 = Ec(1):dE1:Ec(1)+0.1 ;
+        E2 = E1(end):dE2:max(Vtot+DCBO);
+        En=sort([E1 E2]);
+    end
+
+    EEn     = repmat(En' , [1 length(z) ]);
+    V0mat   = repmat(V0  , [length(En) 1]);
+    Vtotmat = repmat(Vtot, [length(En) 1]);
+    Egmat   = repmat(Eg  , [length(En) 1]);
+    meffmat = repmat(meff ,[length(En) 1]);
+    alphamat= repmat(alpha, [length(En) 1]);
+    melin   = meffmat.*(1+alphamat.*(EEn-Vtotmat)); % => linearized mass
+    
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    % Paul Harrisson
+    % Quantum Wells, Wires and Dots.
+    % 4th edition (2016),
+    % chap 2 : "Solutions to Schrodinger's equation"
+    % 2.42: "Two-dimensional systems" page 31
+    
+    meffro = meffmat.*(1+2*alphamat.*EEn); % => linearized mass
+    ro=[];
+    for i=1:length(Ec)
+        ro( En>Ec(i),:,i) = e*meffro(En>Ec(i),:) * m0/(pi*hbar^2);
+        ro( En<Ec(i),:,i) = 0;
+    end
+    
+    % here is what make the code very slow
+    if Ntot==0
+        idx=find(min(CBOt)==CBOt);
+        Ef=-Egt(idx(1))/2;
+        NN=0*Ec';
+        roEf=ro*0;
+    else
+        [Ef,NN,roEf]=find_Ef_f(z,Ec,psic,En,ro,Ntot,T);
+    end
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    
+    ntot2 = repmat(NN,[length(z) 1]).*abs(psic).^2 ;
+    
+    ntot = sum(ntot2,2)' - Dop;  % remove the charge positives (ions)
+        
+    %%%%%%%%%%%%%%%%%%%%%%%%%% Electrical Field %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    F = e*cumtrapz(z,ntot)./(Epsi0*Epsi);
+    MF  = trapz(z,F)/(z(end)-z(1));  % MF=mean(F) on a nonlinear grid z
+    F = F-MF;
+    %%%%%%%%%%%%%%%%%%%%%%%%%%% New Potential %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+    Vsold = Vs;
+    Vs    = -cumsum(F)*dz;
+    %%%%%%%%%%%%%%%%%%%%%%%%%%%% Convergence %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+    Err = abs(  1 - sumVtotVec(end)/sum(Vs)  );
+    sumVtotVec(nloop) = sum(Vs);
+    ErrVec = [ErrVec Err];
+
+    nloop=nloop+1;
+
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Display Results %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+scale_PSI;
+computesISBdipoles;
+
+if DisplayResults == 1
+  PrintResults;
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figures %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+%figure('position',[-3500 100 1000 700],'color','w');
+figure('position',[10 100 1000 700],'color','w');
+subplot(1,1,1,'fontsize',15)
+hold on; grid on; box on;
+col=colormap(jet);
+
+if plot_density==1
+    grid off
+    pcolor(z*1e9,En,sum(ROEf,3)*1e-6 )
+    set(gca,'color',col(1,:))
+    shading flat
+    hcb=colorbar;
+    caxis([0 max(max(sum(ROEf,3)*1e-6))+1])
+    title(hcb,'\fontsize{8}cm-3')
+    
+    if StrainModel == 0
+        plot(z*1e9,V0,  'w--','linewidth',1)
+        plot(z*1e9,Vtot,'w-' ,'linewidth',1)
+    else
+        plot(z*1e9,V0+DCBO,  'k--','linewidth',1)
+        plot(z*1e9,Vtot+DCBO,'k-' ,'linewidth',1)
+    end
+elseif plot_density==0
+    if StrainModel == 0
+        plot(z*1e9,V0,  'b--','linewidth',1)
+        plot(z*1e9,Vtot,'b-' ,'linewidth',1)
+    else
+        plot(z*1e9,V0+DCBO,  'k--','linewidth',1)
+        plot(z*1e9,Vtot+DCBO,'k-' ,'linewidth',1)
+    end
+end
+
+for i=1:length(Ec)
+    plot(z*1e9,PSIc(:,i),'color','r','linewidth',1)
+end
+
+plot([z(1) z(end)]*1e9,[1 1]*Ef,'g','linewidth',1)
+text(z(end)*1e9*0.95,Ef+0.01,'\color{green}Fermi')
+
+xlabel('z (nm)')
+ylabel('Energy (eV)')
+xlim([0 z(end)*1e9])
+if StrainModel == 1
+  title(strcat('T=',num2str(T),'K ; dz=',num2str(dz*1e9),'nm; Ntot=',num2str(Ntot*1e-4,'%.1e'),'cm-2; with STRAIN'))
+else
+  title(strcat('T=',num2str(T),'K ; dz=',num2str(dz*1e9),'nm; Ntot=',num2str(Ntot*1e-4,'%.1e'),'cm-2; without STRAIN'))
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+if plot_convergence==1
+    figure('color','w')
+    semilogy(1:nloop,ErrVec,'bo-')
+    hold on; grid on; box on;
+    xlabel('Cycles')
+    ylabel('Convergence (norm. units)')
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+if plot_field==1
+    
+    figure('color','w')
+    hold on; grid on; box on;
+    [AX,H1,H2]=plotyy(z*1e9,F*1e-2*1e-3,z*1e9,Dop*1e-18*1e-6);
+        
+    set(H1,'color','r')
+    set(H2,'color','b')
+    
+    xlabel('z (nm)')
+    ylabel(AX(1),'E- field (kV/cm)','color','red')
+    ylabel(AX(2),'Doping (1e18 cm-3)','color','blue')
+    
+    set(AX(1),'ycolor','red')
+    set(AX(2),'ycolor','blue')
+        
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+if plot_Vbending==1
+
+    figure('color','w')
+    hold on; grid on; box on;
+
+    [AX,H1,H2]=plotyy(z*1e9,Vbending,z*1e9,ntot*1e-18*1e-6);
+
+    set(H1,'color','r')
+    set(H2,'color','b')
+    
+    xlabel('z (nm)')
+    ylabel(AX(1),'Vbending (eV)','color','red')
+    ylabel(AX(2),'ntot (1e18 cm-3)','color','blue')
+    
+    set(AX(1),'ycolor','red')
+    set(AX(2),'ycolor','blue')
+
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+if plot_mass==1
+    figure('color','w')
+    
+    subplot(1,1,1)
+    hold on; grid on; box on;
+    pcolor(z*1e9,En,melin)
+    
+    colormap(jet)
+    hcb=colorbar;
+    title(hcb,'\fontsize{8}meff')
+    shading flat
+    xlabel('z (nm)')
+    ylabel('Energy (eV)')
+    title('me-lin')
+    xlim([0 z(end)*1e9])
+    
+%     figure
+%     hold on;grid on;
+%     idx=250;
+%     plot(En,melin(:,idx),'r')
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+if plot_ro==1
+    figure('color','w')
+    hold on; grid on; box on;
+    surf(z*1e9,En,sum(ro,3))
+    view(45,30)
+    colormap(jet)
+    hcb=colorbar;
+    title(hcb,'\fontsize{8}ro')
+    shading flat
+    xlabel('z (nm)')
+    ylabel('Energy (eV)')
+    xlim([0 z(end)*1e9])
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+if plot_Epsi==1
+    figure('color','w')
+    hold on; grid on; box on;
+    plot(z*1e9,Epsi,'.-')
+    xlabel('z (nm)')
+    ylabel('Epsilon')
+    xlim([0 z(end)*1e9])
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
diff --git a/TernaryAlloy.m b/TernaryAlloy.m
new file mode 100644
index 0000000..1d0a36a
--- /dev/null
+++ b/TernaryAlloy.m
@@ -0,0 +1,122 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%% III-V Ternary alloys on InP %%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+x=0.53;
+InGaAs = x*InAs + (1-x)*GaAs - x*(1-x)*InGaAs_bowing;   % In0.53Ga0.47As is lattice matched on InP
+x=0.37;
+InGaAs37 = x*InAs + (1-x)*GaAs - x*(1-x)*InGaAs_bowing;   % In0.37Ga0.63As STRAINED on InP
+x=0.45;
+InGaAs45 = x*InAs + (1-x)*GaAs - x*(1-x)*InGaAs_bowing;   % In0.45Ga0.55As STRAINED on InP
+x=0.715;
+InGaAs71 = x*InAs + (1-x)*GaAs - x*(1-x)*InGaAs_bowing;   % In0.715Ga0.285As STRAINED on InP
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+x=0.48;
+AlInAs = x*AlAs + (1-x)*InAs - x*(1-x)*AlInAs_bowing;   % Al0.48In0.52As is lattice matched on InP
+x=0.80;
+AlInAs80 = x*AlAs + (1-x)*InAs - x*(1-x)*AlInAs_bowing;   % Al0.80In0.20As STRAINED on InP
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+x=0.51;
+GaAsSb = x*GaAs + (1-x)*GaSb - x*(1-x)*GaAsSb_bowing;   % GaAs0.51Sb0.49 is lattice matched on InP
+% QCL: Al-fee by G. Strasser's group in Vienna, JVSTB, 32. 02C104 (2014)
+% I did adjust VBO_bowing in order to get CBO=0.36eV between InGaAs/GaAsSb
+
+x=0.56;
+AlAsSb = x*AlAs + (1-x)*AlSb - x*(1-x)*AlAsSb_bowing;   % AlAs0.56Sb0.44 is lattice matched on InP
+% the X-valley is 0.5eV above Gamma point (indirect bandgap material)
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%% III-V Ternary alloys on GaAs %%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+x=0.10;
+AlGaAs_bowing(idx_Eg6c) = -0.127 + 1.310*x;   % from Vurgaftman
+AlGaAs_bowing(idx_VBO)  =  0.27  -  1.05*x;   % to adjust VBO in order to fit CBO from THz QCL
+AlGaAs10 = x*AlAs + (1-x)*GaAs - x*(1-x)*AlGaAs_bowing;
+
+x=0.15;
+AlGaAs_bowing(idx_Eg6c) = -0.127 + 1.310*x;   % from Vurgaftman
+AlGaAs_bowing(idx_VBO)  =  0.27  -  1.05*x;   % to adjust VBO in order to fit CBO from THz QCL
+AlGaAs15 = x*AlAs + (1-x)*GaAs - x*(1-x)*AlGaAs_bowing;
+
+x=0.20;
+AlGaAs_bowing(idx_Eg6c) = -0.127 + 1.310*x;   % from Vurgaftman
+AlGaAs_bowing(idx_VBO)  =  0.27  -  1.05*x;   % to adjust VBO in order to fit CBO from THz QCL
+AlGaAs20 = x*AlAs + (1-x)*GaAs - x*(1-x)*AlGaAs_bowing;
+
+x=0.30;
+AlGaAs_bowing(idx_Eg6c) = -0.127 + 1.310*x;   % from Vurgaftman
+AlGaAs_bowing(idx_VBO)  =  0.27  -  1.05*x;   % to adjust VBO in order to fit CBO from THz QCL
+AlGaAs30 = x*AlAs + (1-x)*GaAs - x*(1-x)*AlGaAs_bowing;
+
+x=0.40;
+AlGaAs_bowing(idx_Eg6c) = -0.127 + 1.310*x;   % from Vurgaftman
+AlGaAs_bowing(idx_VBO)  =  0.27  -  1.05*x;   % to adjust VBO in order to fit CBO from THz QCL
+AlGaAs40 = x*AlAs + (1-x)*GaAs - x*(1-x)*AlGaAs_bowing;
+
+x=0.50;
+AlGaAs_bowing(idx_Eg6c) = -0.127 + 1.310*x;   % from Vurgaftman
+AlGaAs_bowing(idx_VBO)  =  0.27  -  1.05*x;   % to adjust VBO in order to fit CBO from THz QCL
+AlGaAs50 = x*AlAs + (1-x)*GaAs - x*(1-x)*AlGaAs_bowing;
+
+x=0.05;
+InGaAs05 = x*InAs + (1-x)*GaAs - x*(1-x)*InGaAs_bowing;
+x=0.06;
+InGaAs06 = x*InAs + (1-x)*GaAs - x*(1-x)*InGaAs_bowing;
+x=0.07;
+InGaAs07 = x*InAs + (1-x)*GaAs - x*(1-x)*InGaAs_bowing;
+x=0.08;
+InGaAs08 = x*InAs + (1-x)*GaAs - x*(1-x)*InGaAs_bowing;
+x=0.09;
+InGaAs09 = x*InAs + (1-x)*GaAs - x*(1-x)*InGaAs_bowing;
+x=0.10;
+InGaAs10 = x*InAs + (1-x)*GaAs - x*(1-x)*InGaAs_bowing;
+x=0.15;
+InGaAs15 = x*InAs + (1-x)*GaAs - x*(1-x)*InGaAs_bowing;
+x=0.20;
+InGaAs20 = x*InAs + (1-x)*GaAs - x*(1-x)*InGaAs_bowing;
+x=0.25;
+InGaAs25 = x*InAs + (1-x)*GaAs - x*(1-x)*InGaAs_bowing;
+x=0.30;
+InGaAs30 = x*InAs + (1-x)*GaAs - x*(1-x)*InGaAs_bowing;
+x=0.40;
+InGaAs40 = x*InAs + (1-x)*GaAs - x*(1-x)*InGaAs_bowing;
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%% III-V Ternary alloys on InAs/GaSb/AlSb %%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+x=0.25;
+GaInSb25 =  x*InSb + (1-x)*GaSb - x*(1-x)*GaInSb_bowing;
+x=0.35;
+GaInSb35 =  x*InSb + (1-x)*GaSb - x*(1-x)*GaInSb_bowing;
+
+x=0.09;
+GaAsSb91 = x*GaAs + (1-x)*GaSb - x*(1-x)*GaAsSb_bowing;   % GaAs0.09Sb0.91 is lattice matched on InAs
+
+x=0.91;
+InAsSb   = x*InAs + (1-x)*InSb - x.*(1-x)*InAsSb_bowing;    % In0.09As0.91Sb is lattice matched on GaSb, Eg~0.3eV
+
+x=0.5;y=0.08;
+AlGaAsSb = x*GaSb + (1-x)*(y*AlAs + (1-y)*AlSb - y*(1-y)*AlAsSb_bowing); % high barrier lattice matched on GaSb
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%% IV-IV Ternary alloys %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+x=0.25;
+SiGe25 = x*Ge + (1-x)*Si - x*(1-x)*SiGe_bowing;
+x=0.3;
+SiGe30 = x*Ge + (1-x)*Si - x*(1-x)*SiGe_bowing;
+x=0.35;
+SiGe35 = x*Ge + (1-x)*Si - x*(1-x)*SiGe_bowing;
+x=0.4;
+SiGe40 = x*Ge + (1-x)*Si - x*(1-x)*SiGe_bowing;
+x=0.5;
+SiGe50 = x*Ge + (1-x)*Si - x*(1-x)*SiGe_bowing;
+x=0.8;
+SiGe80 = x*Ge + (1-x)*Si - x*(1-x)*SiGe_bowing;
diff --git a/computesISBdipoles.m b/computesISBdipoles.m
new file mode 100644
index 0000000..55c0a60
--- /dev/null
+++ b/computesISBdipoles.m
@@ -0,0 +1,30 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%% Computes ISB dipoles %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% Take care! Some people use meff inside the oscillator strenght f
+% Actually, meff has sens in an infinite QW because there is a single mass value
+% but NOT in multi-QW structure with various materials
+% https://www.nextnano.com/nextnano3/tutorial/1Dtutorial_IntrabandTransitions.htm
+% https://www.nextnano.com/nextnano3/tutorial/1Dtutorial_InGaAs_MQWs.htm
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+z_dipole_c = zeros(length(Ec),length(Ec));
+f_dipole_c = zeros(length(Ec),length(Ec));
+
+for i=1:length(Ec)
+  for j=1:length(Ec)
+    if j>i
+      z_dipole_c(i,j) = abs(  trapz( z , psic(:,i).*z'.*psic(:,j) )  );
+      f_dipole_c(i,j) = 2*m0/hbar^2 * ( Ec(j)-Ec(i) )* e * z_dipole_c(i,j)^2 ;
+    end
+  end
+end
+
+z_dipole_c = z_dipole_c + z_dipole_c.';
+f_dipole_c = f_dipole_c + f_dipole_c.';
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\ No newline at end of file
diff --git a/find_Ef_f.m b/find_Ef_f.m
new file mode 100644
index 0000000..22e8bb4
--- /dev/null
+++ b/find_Ef_f.m
@@ -0,0 +1,54 @@
+function[Ef,NN,roEf]=find_Ef_f(z,Ec,psic,E,ro,Ntot,T)
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Constants %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+e  = 1.602176487E-19;           %% electron charge [C]
+kB = 1.3806488E-23;             %% Boltzmann's constant [J/K]
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+EE   = repmat(E'   ,[1 length(z)]);
+   
+%%%%%%%%%%%%%%%%%%% Computes the Fermi level at any T %%%%%%%%%%%%%%%%%%%%%%%%%%
+if T==0
+   T=1e-10; 
+end
+ 
+Ef=Ec(1);
+Fermi= 1./(1+exp((EE-Ef)/(kB*T/e))) ;
+FFermi=repmat(Fermi,[1 1 length(Ec)]);
+roEf = ro.*FFermi;
+NNz = squeeze(trapz(E,roEf,1));
+NN = trapz(z,NNz.*abs(psic).^2 ,1)  ;
+NtotX=sum(NN);
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%     
+% Now, it will try to get as close as posible to the real Ef with an
+% error of 0.1% by dichotomy
+ddE=0.01; % eV
+Ef1=Ef;
+Ef2=Ef1+ddE;
+
+while  abs(NtotX - Ntot)/Ntot > 0.001  % find the Fermi level at any temperature
+    
+    if NtotX > Ntot
+        Ef  = Ef - abs(Ef1-Ef2)/2 ;  
+        Ef1 = Ef ;
+    else
+        Ef  = Ef + abs(Ef1-Ef2)/2 ;
+        Ef2 = Ef ; 
+    end
+    
+    Fermi  = 1./(1+exp((EE-Ef)/(kB*T/e))) ;  % Fermi Dirac distribution function
+    FFermi = repmat(Fermi,[1 1 length(Ec)]);
+    roEf   = ro.*FFermi;
+    NNz    = squeeze(trapz(E,roEf,1));
+    NN     = trapz(z,NNz.*abs(psic).^2 ,1) ;
+    NtotX  = sum(NN);
+
+end
+
+end
\ No newline at end of file
diff --git a/input_file.m b/input_file.m
new file mode 100644
index 0000000..73aadf1
--- /dev/null
+++ b/input_file.m
@@ -0,0 +1,46 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%% Layers Structure %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% first column is the material used from the "library"
+% second column is the length of the layer in nm = 1e-9m
+% third column is the doping of the layer in 1e18 cm-3 
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% substrate=InP;
+% F0=0;%6e6;%0;               % Electric field [Volt/meter]
+% 
+% M=[
+% AlInAs      5   0
+% AlInAs      1   7
+% AlInAs      5   0
+% InGaAs     6   0
+% AlInAs      5   0
+% InGaAs     7   0
+% AlInAs      5   0
+% AlInAs      1   5
+% AlInAs      5   0
+% ];
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+% substrate=GaAs;      % Important for the Strain model (Si, GaAs, InP, InAs, GaSb)
+% F0=0;%6e6;%0;               % Electric field [Volt/meter]
+% M=[
+% GaAs      10  0.5
+% InGaAs20  15  0
+% GaAs      10  0.5
+% ];
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+substrate=GaAs;      % Important for the Strain model (Si, GaAs, InP, InAs, GaSb)
+F0=0;%16e6;%0;               % Electric field [Volt/meter]
+M=[
+AlGaAs40      5  2
+GaAs         10  0
+AlGaAs40      5  2
+];
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\ No newline at end of file
diff --git a/materialDB_ZB.csv b/materialDB_ZB.csv
new file mode 100644
index 0000000..7423dfd
--- /dev/null
+++ b/materialDB_ZB.csv
@@ -0,0 +1,34 @@
+Material,AlAs,GaAs,InAs,AlSb,GaSb,InSb,AlP,GaP,InP,GaN,AlN,InN,C,Si,Ge,Sn,3C-SiC,ZnSe,ZnTe,CdS,ZnS,CdSe,CdTe,AlGaAs,InGaAs,AlInAs,GaInP,AlInP,AlGaP,GaInSb,AlInSb,AlGaSb,GaAsSb,InAsSb,AlAsSb,GaAsP,InAsP,AlAsP,GaPSb,InPSb,AlPSb,SiGe,Void
+a,5.6611,5.6533,6.0583,6.1355,6.0959,6.4794,5.4672,5.4505,5.8697,4.5,4.38,4.98,3.5668,5.431,5.658,6.49,4.3596,5.668,6.101,5.832,5.409,6.477,6.1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1
+Eg6cG,3.099,1.5192,0.417,2.386,0.812,0.2352,3.63,2.886,1.4236,3.299,5.4,0.78,7.3,4.185,0.89,0.0001,6,2.8215,2.394,2.5,3.68,1.714,1.474,0,0.477,0.7,0.65,-0.48,0,0.415,0.43,0,1.44,0.67,0.8,0.19,0.1,0.22,2.7,1.9,2.7,0,0
+alphaG,8.85E-04,5.41E-04,2.76E-04,4.20E-04,3.00E-04,3.20E-04,5.77E-04,5.77E-04,3.63E-04,5.93E-04,5.93E-04,2.45E-04,0,0.00E+00,5.82E-04,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
+betaG,530,204,93,140,140,170,372,372,162,600,600,624,1,300,296,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1
+EgX,2.24,1.981,1.433,1.696,1.141,0.63,2.52,2.35,2.384,4.52,4.9,2.51,5.5,1.166,1.2,nan,2.36,nan,nan,nan,nan,nan,nan,0.055,1.4,0,0.2,0.38,0.13,0.33,0,0,1.2,0.6,0.28,0.24,0.27,0.22,2.7,1.9,2.7,0,0
+alphaX,7.00E-04,4.60E-04,2.76E-04,3.90E-04,4.75E-04,0.00E+00,3.18E-04,5.77E-04,3.70E-04,5.93E-04,5.93E-04,2.45E-04,0,4.73E-04,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
+betaX,530,204,93,140,94,1,588,372,1,600,600,624,1,636,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1
+EgL,2.46,1.815,1.133,2.329,0.875,0.93,3.57,2.72,2.014,5.59,9.3,5.82,nan,2.4,0.7437,nan,4.6,nan,nan,nan,nan,nan,nan,0,0.33,0,1.03,0,0,0.4,0,0,1.2,0.6,0.28,0.16,0.27,0.22,2.7,1.9,2.7,0,0
+alphaL,6.05E-04,6.05E-04,2.76E-04,5.80E-04,5.97E-04,0.00E+00,3.18E-04,5.77E-04,3.63E-04,5.93E-04,5.93E-04,2.45E-04,0,0,4.77E-04,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
+betaL,204,204,93,140,140,1,588,372,162,600,600,624,1,1,235,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1
+VBO,-1.33,-0.8,-0.59,-0.41,-0.03,0,-0.74,-1.27,-0.986,-2.64,-3.44,-2.34,nan,0,0.638,nan,nan,nan,nan,nan,nan,nan,nan,0,-0.38,-0.64,0,0,0,-0.32,0,0,-0.77,0,-1.71,0,0,0,0,0,0,0,0
+Epsi_stat,10.6,12.9,15.15,12.04,15.7,16.8,9.8,11.1,12.5,9.7,9,15,nan,11.7,16,nan,nan,9.2,9.3,8.6,8.6,9.5,2.27,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
+Epsi_inf,8.16,10.6,12.3,12.3,14.5,15.7,nan,9.11,10.9,5.3,4.8,8.5,nan,nan,nan,nan,nan,5.8,6.9,5.3,5.2,6.2,nan,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
+EP_Kane,21.1,22.71,21.11,18.7,22.88,22.49,nan,nan,17,nan,nan,nan,0,0,0,0,0,nan,nan,nan,nan,nan,nan,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
+EP_Luttinger,21.1,28.8,21.5,18.7,27,23.3,17.7,31.4,20.7,25,27.1,17.2,0,0,0,0,0,24.2,19.1,21,20.4,20,20.7,0,-1.48,-4.81,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
+Dso,0.28,0.341,0.39,0.676,0.752,0.81,0.07,0.08,0.108,0.017,0.019,0.005,0.006,0.044,0.297,0,0.01,0.403,0.97,0.08,0.089,0.441,0.949,0,0.15,0.15,0,-0.19,0,0.1,0.25,0.3,0.6,1.2,0.15,0,0.16,0,0,0.75,0,0,0
+F,-0.48,-1.94,-2.9,-0.56,-1.63,-0.23,-0.65,-1.97,-1.31,-0.95,-1.01,-4.36,-1,-1,-1.54,-1,-1,-0.76,-0.007,-1.08,-1.78,-1.39,-1.55,0,1.77,-4.44,0.78,0,0,-6.84,0,0,0,0,0,0,0,0,0,0,0,0,0
+gamma1,3.76,6.98,19.49,5.18,13.4,34.8,3.35,4.05,5.08,2.7,1.92,3.72,4.24,4.26,13.26,-15,2.817,3.77,3.74,4.11,2.54,5.51,5.2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
+gamma2,0.82,2.06,8.31,1.19,4.7,15.5,0.71,0.49,1.6,0.76,0.47,1.26,0.82,0.38,4.24,-11.5,0.508,1.24,1.07,0.77,0.75,1.24,1.89,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
+gamma3,1.42,2.93,9.14,1.97,6,16.5,1.23,2.93,2.1,1.11,0.85,1.63,1.71,1.56,5.64,-8.6,0.86,1.67,1.64,1.53,1.09,2.14,2.46,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
+Eg8c,4.7,4.66,4.577,3.77,3.69,3.7752,4.82,4.47,4.97,nan,nan,nan,nan,3.4,3.31,nan,nan,nan,nan,nan,nan,nan,nan,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
+Eg7c,4.55,4.49,4.387,3.53,3.35,3.3852,4.78,4.38,4.78,nan,nan,nan,nan,3.4,3.12,nan,nan,nan,nan,nan,nan,nan,nan,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
+Eg6v,-11.73,-12.5,-12.64,-10.62,-11.47,-10.54,-11.21,-12.14,-11.04,nan,nan,nan,nan,-12.5,-12.64,nan,nan,nan,nan,nan,nan,nan,nan,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
+ac,-5.64,-7.17,-5.08,-4.5,-7.5,-6.94,-5.7,-8.2,-6,-6.71,-4.5,-2.65,nan,-5.1,-9.5,nan,nan,nan,nan,nan,nan,nan,nan,0,2.61,-1.4,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
+av,-2.47,-1.16,-1,-1.4,-0.8,-0.36,-3,-1.7,-0.6,-0.69,-4.9,-0.7,nan,0,0,nan,nan,nan,nan,nan,nan,nan,nan,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
+bv,-2.3,-2,-1.8,-1.35,-2,-2,-1.5,-1.6,-2,-2,-1.7,-1.2,nan,2.33,2.08,nan,nan,nan,nan,nan,nan,nan,nan,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
+dv,-3.4,-4.8,-3.6,-4.3,-4.7,-4.7,-4.6,-4.6,-5,-3.7,-5.5,-9.3,nan,0,0,nan,nan,nan,nan,nan,nan,nan,nan,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
+c11,1250,1221,832.9,876.9,884.2,684.7,1330,1405,1011,293,304,187,nan,1675,1315,nan,nan,nan,nan,nan,nan,nan,nan,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
+c12,534,566,452.6,434.1,402.6,373.5,630,620.3,561,159,160,125,nan,650,494,nan,nan,nan,nan,nan,nan,nan,nan,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
+c44,542,600,395.9,407.6,432.2,311.1,615,703.3,456,155,193,86,nan,796,0,nan,nan,nan,nan,nan,nan,nan,nan,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
+me,0.15,0.067,0.023,0.14,0.042,0.014,0.22,0.79,0.079,0.193,0.316,0.054,nan,1.08,0.22,nan,nan,0.13,0.122,0.14,0.34,0.11,0.099,0,0.0091,0.032,0.051,0.22,0,0.0092,0,0,0.049,0.035,0,0,0,0,0,0,0,0,0
+mhh,0.76,0.49,0.41,0.5,0.3,0.4,0.7,0.79,0.6,0.81,1.33,0.835,nan,0.81,0.34,nan,nan,nan,nan,nan,nan,nan,nan,0,-0.145,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0
+mlh,0.16,0.08,0.026,0.11,0.05,0.02,nan,0.14,0.12,0.265,0.466,0.08,nan,nan,nan,nan,nan,nan,nan,nan,nan,nan,nan,0,0.0202,0,0,0,0,0.011,0,0,0,0,0,0,0,0,0,0,0,0,0
diff --git a/scale_PSI.m b/scale_PSI.m
new file mode 100644
index 0000000..582b409
--- /dev/null
+++ b/scale_PSI.m
@@ -0,0 +1,17 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%% Scaling and shifting the wavefunctions %%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+for i=1:length(Ec)
+  PSIc(:,i)=abs(psic(:,i)).^2/max(abs(psic(:,i)).^2)*ScF + Ec(i); % normalisation for the plotting
+end
+
+% ROEf is only build to plot the density of electron in a 2D map
+for ii=1:length(Ec)
+    psic_M = repmat(psic(:,ii),[1,length(En)])';
+    ROEf(:,:,ii) = roEf(:,:,ii) .* abs(psic_M.^2);
+end
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\ No newline at end of file