NaN error when i add P =P + delta_P:

% Constants
L=10;
Nx = 6; % Number of grid points
dx = L / (Nx - 1); % Grid spacing  
x = linspace(0, L, Nx); % Spatial points
% Constants from your provided values
f = 37.40e-3;
c = 343;
D = 104e-3;
g = 9.81;
epsilon = 1e-6;
initial_guess_m = 57.88;  % Initial guess for m
initial_P_guess = 1.015e5; % Initial guess for P

% Initialize arrays for m and P with initial guesses
m = initial_guess_m*ones(Nx, 1);
P = initial_P_guess*ones(Nx, 1);

% Set boundary conditions
P(1) = 1.13e6; % Boundary condition at i=1
boundary_condition_mNx = 57.88;  % Boundary condition at m(Nx)

% Convergence tolerance
tolerance = 1e-6;

% Maximum number of iterations
max_iterations = 100;

% Initialize variables
iteration = 0;
converged = false;
m(Nx)=57.88;

% Newton-Raphson iteration
while ~converged && iteration < max_iterations
 
    % Section - Residual calculation of continuity & Momentum
    Rm = continuity_residual(m, Nx, dx, boundary_condition_mNx);
    RP = momentum_residual(P, m, dx, f, c, D,Nx);
  
    % Section - Jacobian Matrix calculation for continuity and momentum
    Jmm =continuity_jacobian(m, Nx, dx);
    JP = momentum_jacobian(P, m, dx, f, c, D, Nx);

    % Section - Combined Residual
    combined_residual = [Rm; RP];

    % Section - Combined Jacobian Matrix
    combined_Jacobian = [Jmm,zeros(Nx, Nx);zeros(Nx, Nx),JP];
    
    % Section - Solve Linear system
    delta_combined = - combined_Jacobian \ combined_residual;
   
    % Section - Split the combined delta into delta_P and delta_m
    delta_m = delta_combined(1:Nx);
    delta_P = delta_combined(Nx+1:end);
   
    % Update the variables
    m = m + delta_m ;
    P =P + delta_P;
    
    % Check for convergence
    if max(abs(delta_combined )) < tolerance
        converged = true;
    end
    
    % Increment the iteration counter
    iteration = iteration + 1;
end

% Plot pressure m
figure;
plot(x, m, 'r');  % 
xlabel('Spatial Position (x)');
ylabel('mass flowrate (kg/s)');
title('mass flowrate  Distribution');

% Plot pressure P
figure;
plot(x, P, 'b');  % 
xlabel('Spatial Position (x)');
ylabel('Pressure (Pasc)');
title('Pressure Distribution');

%% Residual vector 
% Function for continuity equation residual
function Rm = continuity_residual(m, Nx, dx, boundary_condition_mNx)
    Rm = zeros(Nx, 1);
    
    % Calculate the residual for interior points
    for i = 2:Nx - 1
        Rm(i) = (m(i + 1) - m(i - 1)) / (2 * dx);%% Central discretization
    end
    
    % Handle residual at i = 1
    Rm(1) = (m(2) - m(1)) / dx; %% Forward discretization to avoid ghost point i=0
    
    % Handle boundary condition at i = Nx
    Rm(Nx) = m(Nx)-boundary_condition_mNx; 
end

% Function for momentum equation residual
function RP = momentum_residual(P, m, dx, f, c, D, Nx)
    % Calculate residuals for the interior grid points
    RP = zeros(Nx, 1);
     
    % Handle boundary condition at i = 1
    RP(1) = P(1)-1.13e6;
    for i = 2:Nx - 1
        RP(i) = (P(i+1) - P(i-1)) / (dx) + f * (m(i)^2 * c^2 ) /2*( P(i)) ;
    end
     % Include Residual momentum at i=Nx (Backward discrteization)
     RP(Nx) = (P(Nx) - P(Nx-1)) / (dx) +  (f *m(Nx)^2 * c^2) / 2*(P(Nx));% Backward discretization to avoid ghost point i=Nx+1
         
end


%% Jacbian Matrix 
% Continuity Jacobian Function
function Jmm = continuity_jacobian(m, Nx, dx)
    Jmm = zeros(Nx, Nx);
    
    % Calculate the Jacobian matrix for interior points
    for i = 2:Nx - 1
        Jmm(i, i - 1) = -1 / (2 * dx);
        Jmm(i, i + 1) = 1 / (2 * dx);
    end
    Jmm(1,1) = -1/dx;%derivative of dRm(1)/dm(1)
    Jmm(1,2) = 1/dx;%derivative of dRm(1)/dm(2)
    Jmm(Nx,Nx) = 1.0;%derivative of dRm(Nx)/dm(Nx) at the boundary condition point m(Nx)
end

function JP = momentum_jacobian(P, m, dx, f, c, D, Nx)
    % Initialize the Jacobian matrix
    JP = zeros(Nx, Nx);

    % Define the main diagonal and upper/lower diagonals
    main_diag = zeros(Nx, 1);
    upper_diag = zeros(Nx-1, 1);
    lower_diag = zeros(Nx-1, 1);
    
    % Calculate values for the diagonals
    main_diag(2:Nx-1) = - (f * m(2:Nx-1).^2 * c^2 ) ./ (2 * (P(2:Nx-1).^2));
    upper_diag(2:Nx-1) = 1 / (2 * dx);
    lower_diag(2:Nx-1) = -1 / (2 * dx);
  
    % Set the diagonals in the Jacobian matrix
     JP = diag(main_diag) + diag(upper_diag, 1) + diag(lower_diag, -1);
     JP(1, 1) = 1  ;%derivative of dRP(1)/dP(1) at the boundary condition point P(1)
     JP(Nx, Nx-1) = -1/dx;%derivative of dRP(Nx)/dP(Nx-1) at the last grid point Nx with respect to the second-to-last grid point Nx-1
     JP(Nx, Nx) = 1/dx - (f * m(Nx)^2 * c.^2) / (2 * P(Nx).^2);% derivative of dRP(Nx)/dP(Nx) at the last grid point Nx
end

Error : When i add P =P + delta_P; to update the variable my jacobian matrix and residual vector tunrs to NaN error?
    

[Torsten]

Since you access m(i) in the computation of the momentum residuals, the bold part of the combined Jacobian below

combined_Jacobian = [Jmm,zeros(Nx, Nx);zeros(Nx, Nx),JP];

must be wrong.

[MAYUR MARUTI DHIRDE]

If i am running my code with this combined_Jacobian = [Jmm,zeros(Nx, Nx);JP];

Then i get an error ,

Error

Dimensions of arrays being concatenated are not consistent.

Error in (line 51)

combined_Jacobian = [Jmm,zeros(Nx, Nx);JP];

[Torsten]

If i am running my code with this combined_Jacobian = [Jmm,zeros(Nx, Nx);JP];

Sure. I didn't say you should remove this part, but change zeros(Nx,Nx) in order to account for the dependencies of the momentum residuals on the continuity variables m(i).

And I'm not sure since I don't remember your mathematical equations, but I think it should be

for i = 2:Nx - 1
    RP(i) = (P(i+1) - P(i-1)) / (dx) + f * (m(i)^2 * c^2 ) /(2* P(i)) ;
end
% Include Residual momentum at i=Nx (Backward discrteization)
RP(Nx) = (P(Nx) - P(Nx-1)) / (dx) +  (f *m(Nx)^2 * c^2) / (2*P(Nx));% Backward discretization to avoid ghost point i=Nx+1

instead of

for i = 2:Nx - 1
    RP(i) = (P(i+1) - P(i-1)) / (dx) + f * (m(i)^2 * c^2 ) /2*( P(i)) ;
end
% Include Residual momentum at i=Nx (Backward discrteization)
RP(Nx) = (P(Nx) - P(Nx-1)) / (dx) +  (f *m(Nx)^2 * c^2) / 2*(P(Nx));% Backward discretization to avoid ghost point i=Nx+1