Implicit finite difference for 5 simultaneous pde and Newton-Raphson method

Hello,

I am trying to solve my system with 5 nonlinear pde with 5 unknown functions using implicit finite difference method. At the same time, the code uses Newton-Raphson iteration for gap1_w+gap2_w=1. I have coded the problem as shown below

%--------------------------------------------------------------------------
Nt = 3;                                % Total time step number
Nx = 101;                              % Total grid
dt = 0.0001;                           % Time step
dx = 1/(Nx-1);                         % Grid size
xplot = 1:dx:2;                        % Record the x scale for plots
r = dt/dx;                             % define the ratio r
iplot = 1;                             % Counter used to count plots
tplot(1) = 0;                          % Initial time in plots
Costep = 300;                          % Maximum number of iterations
nplots = 50;                           % Number of snapshots (plots) to take
plot_step = Costep/nplots;             % Number of time steps between plots
dgap = 0.001;                          % Newton-Raphson interval

%--------------------------------------------------------------------------
%Preallocation and initial conditions
gap1_w =.5*ones(1,Nx);    gap2_w =.5*ones(1,Nx);
u1_w = ones(1,Nx);        u2_w = ones(1,Nx);
p1_w = zeros(1,Nx);       gap1_wb = gap1_w;
%Setting up auxillary variables
for i=1:Nx
    im(i)=i-1;
end
im(1)=1;
%--------------------------------------------------------------------------
for nt=1:Nt+1
    t=(nt-1)*dt; 
    gap1_w(1) = .5*(1+.5*sin(pi*t));              % LHS BC for gap1_w
    gap2_w(1) = .5*(1-.5*sin(pi*t));              % LHS BC for gap2_w
    % From the mass conservation H1.u1 +H2.u2 =1 --------------------------
    beta = (2+ cos(pi*t));
    delta = (2+ cos(pi*t));
    B = (gap1_w(1)*beta)+(gap2_w(1)*delta);
    epsilon = 1/B;
    %----------------------------------------------------------------------
    u1_w(1) = epsilon*(2+ cos(pi*t));             % LHS BC for u1_w
    u2_w(1) = epsilon*(2+ cos(pi*t));             % LHS BC for u2_w
    p1_w(1) = epsilon*pi*(sin(pi*t));             % LHS BC for p1_w
    %----------------------------------------------------------------------
    % Save the values from the previous time-step as they are needed for the time derivative finite difference- 
    gap1n=gap1_w;
    gap2n=gap2_w;
    u1n=u1_w;
    u2n=u2_w;
    p1n=p1_w;
    % Newton - Raphson for gap1_w + gap2_w = 1-----------------------------
    flag=0;
    mmm=1;
    m=1;
    gap1_w = gap1_wb;
    gap2_w = 1-gap1_w;
    while (flag == 0)
        gap2_w = 1-gap1_w;
        % Compute new gap1_w, gap2_w, u1_w, u2_w, p1_w using first order backward time finite difference scheme.
        for i=1:Nx
            % For gap 1----------------------------------------------------
            gap1_w(i) = gap1n(i) - r*(u1_w(i).*(gap1_w(i)-gap1_w(im(i)))+ ...
                        gap1_w(i).*(u1_w(i)-u1_w(im(i))));

              u1_w(i)=u1n(i)-r*(u1_w(i).*(u1_w(i)-u1_w(im(i)))+ (p1_w(i) - p1_w(im(i))));
              % For gap 2----------------------------------------------------
              gap2_w(i) = gap2n(i) - r*(u2_w(i).*(gap2_w(i)-gap2_w(im(i)))+ ...
                          gap2_w(i).*(u2_w(i)-u2_w(im(i))));

              u2_w(i)=u2n(i)-r*(u2_w(i).*(u2_w(i)-u2_w(im(i))+(p1_w(i) - p1_w(im(i)))));
              % Above, we have used p1 instead of p2 as they are equal.
          end
          for i=1:Nx
          st(m,i) = 1-(gap1_w(i)+gap2_w(i));
          end
          m = m + 1;
          if (m == 2)
              gap1_w  = gap1_w + dgap;
              gap2_w = 1- gap1_w;
              st(m)           
          end
          if (m == 3)
              for i=1:Nx
              gap1_w(i)  = (gap1_w(i).*st(1,i)-(gap1_w(i)-dgap).*st(2,i))./(st(1,i)-st(2,i));
              end
              gap2_w = 1- gap1_w;
          end
          if (m == 4)
              mmm = mmm + 1;
              if (mmm == 2)
                  m = 1;
              else
                  flag = 1;
              end
          end
      end
      % Periodically record gap1_w, gap2_w, u1_w, u2_w, p1_w for ploting.

      gap1_wplot(:,iplot) = gap1_w(:);  % Record gap1_w(i) for plotting
      gap2_wplot(:,iplot) = gap2_w(:);  % Record gap2_w(i) for plotting
      total(:,iplot) = gap1_w(:) + gap2_w(:); % Total gap for plotting
      mass(:,iplot) = gap1_w(:).*u1_w(:) + gap2_w(:).*u2_w(:); % Total gap for plotting
      u1_wplot(:,iplot) = u1_w(:);      % Record u1_w(i) for plotting
      u2_wplot(:,iplot) = u2_w(:);      % Record u2_w(i) for plotting
      p1_wplot(:,iplot) = p1_w(:);      % Record p1_w(i) for plotting
      tplot(iplot) = nt*dt;             % Record time for plots
      %         end
      iplot = iplot+1;

end

The thing is, it should use gap1_w,gap2_w u1_w,u2_w, p1_w on the right hand side in these lines using CURRENT TIME STEP NOT PREVIOUS

gap1_w(i) = gap1n(i) - r*(u1_w(i).*(gap1_w(i)-gap1_w(im(i)))+ ...
                        gap1_w(i).*(u1_w(i)-u1_w(im(i))));
              u1_w(i)=u1n(i)-r*(u1_w(i).*(u1_w(i)-u1_w(im(i)))+ (p1_w(i) - p1_w(im(i))));
              % For gap 2----------------------------------------------------
              gap2_w(i) = gap2n(i) - r*(u2_w(i).*(gap2_w(i)-gap2_w(im(i)))+ ...
                          gap2_w(i).*(u2_w(i)-u2_w(im(i))));
              u2_w(i)=u2n(i)-r*(u2_w(i).*(u2_w(i)-u2_w(im(i))+(p1_w(i) - p1_w(im(i)))

However, I assume it uses these from the previous time step. For example, instead of the gap1_w from the CURRENT time step, it recognises gap1_w as gap1n since it did not complete the loop.

Any help is much appreciated!