function [r1,force,air,summ,pres,q,flow_sum]=coup_DM(Tmax,n,order,kappa)
% to reduce to Austin's model, k1,k2,k7=0
% Force= dt and Rref=0.2792

delta_t=Tmax/n;


numBr = 2^order - 1;     % number of airways
rimax = [0.296 0.318 0.337 0.358 0.384 0.414 0.445 0.484 0.539 0.608 ...
    0.692 0.793 0.913 1.052 1.203 1.374];
rinitial=[];
for k=1:order
    rinitial=[rinitial;rimax(k)*ones(2.^(order-k),1)];
end
rinitial=rinitial+0.01*rand(1,length(rinitial))';

r1 = zeros(length(rinitial),n+1);
epsilon = 0;
%% Distribution Moment ICs
% M10=0.005;
% M11=0.001;
% M12=0.01;
% M20=0.005;
% M21=0.001;
% M22=0.01;
% C=1;
% 
%  
  M10= 0.6209;
  M11=  0.3106;
  M12= 0.2071;
  M20= 0.1749;
  M21= 0.0871;
  M22= 0.0579;
  C= 0.222222222222222;
%    
   
RR = [M10;M11;M12;M20;M21;M22;C];
r1(:,1)=rinitial;

%% Initialising the structure array
output.rdot = NaN;
output.R = NaN;
output.Ptm = NaN;
output.tau = NaN;
output.mu = NaN;
output.q = NaN;
output.p = NaN;
output.pbot = NaN;
output.deltap = NaN;
output.W = NaN;
output.lambda = NaN;
output.Lambda = NaN;
output.gammajplus = NaN;
output.gammajminus = NaN;

for i=1:n
%      if mod(i,n) == 100
%         fprintf('\niteration #: %g \n',i)
%     end 
% t(i)=i*delta_t;
%     dt(i)=0.4*t(i);
    parms = [kappa;epsilon];


    [k1,F1,~,~]=Type1p2_Connect_GEN_v3d0_DM(RR,r1(:,i),order,parms);
    [k2,F2,~,~]=Type1p2_Connect_GEN_v3d0_DM(RR+(delta_t/2)*F1,r1(:,i)+(delta_t/2)*k1,order,parms);
    [k3,F3,~,~]=Type1p2_Connect_GEN_v3d0_DM(RR+(delta_t/2)*F2,r1(:,i)+(delta_t/2)*k2,order,parms);
    [k4,F4,out,force(:,i),pbot(i),q(:,i)]=Type1p2_Connect_GEN_v3d0_DM(RR+(delta_t)*F3,r1(:,i)+(delta_t)*k3,order,parms);
    
    
    Rnew = RR + (delta_t/6)*(F1+(2*F2)+(2*F3)+F4);
    RR=Rnew;
    
     r1(:,i+1) = r1(:,i) +(delta_t/6)*(k1+(2*k2)+(2*k3)+k4);

end


time = 0:delta_t:Tmax;
CM=jet(length(rinitial));
air=zeros(1,numBr);
for i=1:length(rinitial)
      hold on
%     subplot(8,4,i)
%   
    plot(time,r1(i,:));
% flow_mean(i)=mean(flow(i,:));
  if r1(i,end) <= 0.01
       air(i)=0;
   else
       air(i)=1;
  end

  
end


numOrd1 = (numBr + 1)/2;
x = cell(1,numOrd1);
y=cell(1,numOrd1);
 for k=1:length(y)
 y{k}=q(k,end);
 end
 maxx=(max([y{:}]));
 % % To normalise the flows so that it is between [0,1]
  for l = 1:length(y)
      y{l} = y{l}./maxx;
  end
%   for j=1:length(y)
%       yy{j}= y{j}./max([y{:}]);
%   end
  flow_sum = zeros(1,length(y{1}));
  for i = 1:length(y)
      switch i
          case 1
              neighbour_diff1 = sqrt((y{i} - y{end}).^2 + (y{i} - y{i+1}).^2);
          case length(y)
              neighbour_diff1 = sqrt((y{i} - y{i-1}).^2 + (y{i} - y{1}).^2);
          otherwise
              neighbour_diff1 = sqrt((y{i} - y{i-1}).^2 + (y{i} - y{i+1}).^2);
      end
     flow_sum = flow_sum + neighbour_diff1;
  end
flow_sum=flow_sum./numOrd1;
%  flow_sum=flow_sum./(2^(order-3)*4*sqrt(2));


for i = 1:length(x)
    x{i} = air(i);
end
summ = zeros(1,length(x{1}));
for i = 1:length(x)
    switch i
        case 1
            neighbour_diff = sqrt((x{i} - x{end}).^2 + (x{i} - x{i+1}).^2);
        case length(x)
            neighbour_diff = sqrt((x{i} - x{i-1}).^2 + (x{i} - x{1}).^2);
        otherwise
            neighbour_diff = sqrt((x{i} - x{i-1}).^2 + (x{i} - x{i+1}).^2);
    end
    summ = summ + neighbour_diff;
    
end
summ=summ./(2^(order-3)*4*sqrt(2));
pres=pbot(end);


