% File: pde_3.m
function yt=programme_BTE(~,y)  
%
% Global area
global nr  ne  dr N  M c1 c2 c6  c  ncall
%
% 1D to 2D matrices
  for i=1:nr
  for j=1:ne
    ij=(i-1)*ne+j;
    N(i,j)=y(ij);    
    
  end
  end
  
  % Preallocation for speed
  frd=zeros(nr,ne);
  ftbd=zeros(nr,ne);
  N1=zeros(nr,ne);
  N2=zeros(nr,ne);
  f=zeros(nr,ne);
 
  Nes1=zeros(nr,ne);
  Nes2=zeros(nr,ne);
  Nes3=zeros(nr,ne);
  Nes4=zeros(nr,ne);
  fes=zeros(nr,ne);
  Nr=zeros(nr,ne);
  Ne=zeros(nr,ne);
  Nt=zeros(nr,ne);
  yt=zeros((nr*ne),1);
  yt=yt';
  
 % Exciton phonon scattering term

for j=(ceil(M/4)+1):ne   
sum=0;
 p1= j-1;
 
if j<=M+1
     p4=1;
 else
p4=((sqrt(j-1)-sqrt(M)).^2)+1;
 end
 p2= ceil(p4);
 
for p=p2:p1
 
sum= sum+((j-p).^2).*((N(:,j).*(1+(1/(exp(c.*(j-p))-1))).*(1+N(:,p)))-((1+N(:,j)).*(1/(exp(c.*(j-p))-1)).*N(:,p)));
end
 
N1(:,j)=(-(1.844)/((M^2)*sqrt(M)*sqrt(j-1))).*sum;
end

for j=2:ne
sum=0;

p3= ((sqrt(j)-sqrt(M)).^2);
p4=((sqrt(j)+sqrt(M)).^2);

p1=max(ceil(p3),(j+1));

p2=min(ne, ceil(p4));
for p=p1:p2
    
  
sum= sum+((j-p).^2).*((N(:,j).*(1/((exp(c*(p-j)))-1)).*(1+N(:,p)))-((1+N(:,j)).*(1+(1/((exp(c*(p-j)))-1))).*N(:,p)));
end
N2(:,j)=((-(1.844)/((M^2)*sqrt(M)*sqrt(j-1))).*sum);        

end
N2(:,1)=((7.376)*((N(:,(M+1)).*(1+(1/((exp(c.*M))-1))))-(N(:,1).*((1/((exp(c.*M))-1))-N(:,(M+1)))))/(1));      


  for i=1:nr
      for j=1:ne
          f(i,j)=(N1(i,j)+N2(i,j));
          
      end
  end
  
 


% Exciton exciton scattering
for j=2:ne
    sum1=0;
   
    for j1=1:j
        for j2=(j-j1):ne
            j3=j1+j2-j;
            if j3>ne
               j3=ne;
            elseif j3<1
               sum1=sum1;
            else 
          
            sum1=sum1+(((N(:,j1).*N(:,j2).*(1+N(:,j)).*(1+N(:,j3)))-(N(:,j).*N(:,j3).*(1+N(:,j1)).*(1+N(:,j2)))).*(min([sqrt(j1-1),sqrt(j2-1),sqrt(j-1),sqrt(j3-1)])));
            end
        end
    end
    Nes1(:,j)=c6*(1/((M^(2))*sqrt(j-1))).*sum1;
end



for j=2:ne
    sum2=0;
   
    for j1=j:ne
        for j2=1:ne
            j3=j1+j2-j;
            if j3>ne
               j3=ne;
            elseif j3<1
               sum2=sum2;
            else 
          
            sum2=sum2+(((N(:,j1).*N(:,j2).*(1+N(:,j)).*(1+N(:,j3)))-(N(:,j).*N(:,j3).*(1+N(:,j1)).*(1+N(:,j2)))).*(min([sqrt(j1-1),sqrt(j2-1),sqrt(j-1),sqrt(j3-1)])));
            end
        end
    end
    Nes2(:,j)=c6*(1/((M^(2))*sqrt(j-1))).*sum2;
end

%at k=0
sum2=0;
   
    for j1=1:ne
        for j2=1:ne
            j3=j1+j2-1;
            if j3>ne
               j3=ne;
            elseif j3<1
               sum2=sum2;
            else 
          
            sum2=sum2+(((N(:,j1).*N(:,j2).*(1+N(:,1)).*(1+N(:,j3)))-(N(:,1).*N(:,j3).*(1+N(:,j1)).*(1+N(:,j2)))));
            end
        end
    end
Nes2(:,1)=c6*(1/((M^(2)))).*sum2;

c7=2.3*(10^(-5));    %[unit=1/ns][c7=(4*((al)^2)*(sqrt(e0)))/(pi*V1*(sqrt(2*Mx)))][a1=scattering length=2 nm, e0=0.6 meV, V1=volume=1 micro_meter^3, Mx=exciton mass=2.61*9.109*10^(-31)kg]  
%at k2=j3=0
for j=1:ne
    sum3=0;
     for j1=1:j
        
            j2=j-j1;
            if j2>ne
               j2=ne;
            elseif j2<1
               sum3=sum3;
            else 
          
            sum3=sum3+((N(:,j1).*N(:,j2))-(N(:,j).*(1+N(:,j1)+N(:,j2))));
            end
       
     end
     Nes3(:,j)=c7*(1/(sqrt(M)*sqrt(j))).*N(:,1).*sum3;            
end
%at p2=j2=0
for j=1:ne
    sum4=0;
     for j1=j:ne
        
            j3=j1-j;
            if j3>ne
               j3=ne;
            elseif j3<1
               sum4=sum4;
            else 
          
            sum4=sum4+((N(:,j1).*(1+N(:,j3)+N(:,j)))-(N(:,j).*N(:,j3)));
            end
       
     end
     Nes4(:,j)=2*c7*(1/(sqrt(M)*sqrt(j))).*N(:,1).*sum4;        
end

   
   
for i=1:nr
    for j=1:ne
        fes(i,j)=Nes1(i,j)+Nes2(i,j)+Nes3(i,j)+Nes4(i,j);
       
    end
end

% radiative decay
for i=1:nr
      for j=1:ne
          frd(i,j)=(1/650)*N(i,j);
      end
end
 

% Two body decay
App=10^(-6);     %[Auger constant, unit=(micro_m^3)/ns] [10^(-6) means App=10^(-18 and 10^(-4) means App=10^(-16))
Cb=52.284*(10^(6));   %[(8*pi*(Mx^(3/2))*(e0^(3/2)))/(sqrt(2)^h1^3)][unit=1/(micro_m)^3]
 i=1:nr;
    sum=0;
    for j=2:ne
        sum=sum+(sqrt(j-1).*N(i,j));
    end
    N_e_sum(i)=(Cb/(M^(3/2)))*sum;


 
    sum1=0;
    for j=2:ne
        sum1=sum1+(Cb/(M^(3/2)))*(sqrt(j-1));
    end



for i=1:nr
      for j=2:ne
          ftbd(i,j)=0.004*(-(App*N_e_sum(i).*N(i,j))+((1/2)*((App))*((N_e_sum(i)).^2).*(1/sum1)));
      end
 end


% Boundary conditions
  for i=1:nr
  for j=1:ne
      
      
      if(i==1)
         Nr(i,j)=0;
      elseif(i==nr)
         Nr(i,j)=0;
      %elseif(j==1)
       %  Nr(i,j)=0; 
      elseif(j==ne)
         Nr(i,j)=0;
      else 
          Nr(i,j)=((N(i,j)-N(i-1,j))/dr);
     end
  
      
      if(j==ne)
          Ne(i,j)=0;
      elseif(j==1)
          Ne(i,j)=0;
      %elseif(i==1)
       %  Ne(i,j)=0;
      elseif(i==nr)
          Ne(i,j)=0;
      else 
          Ne(i,j)=((N(i,j+1)-N(i,j)));
      end
      
      
      
          Nt(i,j)=-(c1*sqrt(j).*Nr(i,j))+(c2*sqrt(j).*(i-1).*dr.*Ne(i,j))+f(i,j)+fes(i,j)-frd(i,j)+ftbd(i,j);
         
  end
  end
  %
% 2D to 1D matrices
  for i=1:nr
  for j=1:ne
    ij=(i-1)*ne+j;
    yt(ij)=Nt(i,j);
  end
  end 
%
% Transpose and count
  yt=yt';
  ncall=ncall+1;