function [r,air,summ,flow_sum,pres,out]=coup(Tmax,n,order,kappa)%,rinitial)
tic
epsilon = 0;

delta_t=Tmax/n;

M=2000;
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

r = zeros(length(rinitial),n+1);
xrange=[-30 30];
X0=linspace(xrange(1),xrange(2),M);

r(:,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 j=1:30
for i=1:n
    
    parms = [kappa;epsilon];

    
     t(i)=i*delta_t;
    [k1,~]=Type1p2_Connect_GEN_v3d0(r(:,i),order,parms);
    [k2,~]=Type1p2_Connect_GEN_v3d0(r(:,i)+(delta_t/2)*k1,order,parms);
    [k3,~]=Type1p2_Connect_GEN_v3d0(r(:,i)+(delta_t/2)*k2,order,parms);
    [k4,out,pbot(i),p(:,i),q(:,i),mu(:,i)]=Type1p2_Connect_GEN_v3d0(r(:,i)+(delta_t)*k3,order,parms);
        
     r(:,i+1) = r(:,i) +(delta_t/6)*(k1+(2*k2)+(2*k3)+k4);
end

time = 0:delta_t:Tmax;
CM=jet(length(rinitial));

for i=1:length(rinitial)

   if r(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;


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);
toc
% profile viewer
%    psi=air(1)+air(2);
%    psi1=air(3)+air(4);
   
%    A_diff=diff([air(1:4) ;air(1)]');
% A_sq=A_diff.^2;
% A_norm1=sqrt(sum(A_sq))./4;
% psi=air(1)+air(2);
%    psi1=air(3)+air(4);
% end
% nn=1:30;
% air_state=[rinitial air'];
for i = 1:length(rinitial)
    plot(time,r(i,:),'-','color',CM(i,:),'LineWidth', 1.5, 'MarkerSize', 8);
    hold on
end
%     
%     
legendStr = cell(1,length(rinitial));
for i = 1:length(legendStr)
    legendStr{i} = sprintf('r_{%d}',i);
%     text(max(time),max(r(i)),'r_{%d}')
% text(10,max(r(i,:)),legendStr{i})
end
legend(legendStr,'Location','BestOutside');
% 
% 
% xlabel('Time (s)');
% ylabel('radius (mm)');
% grid on;