Ode15s , Plotting issue.

 
Nx = 50;         % Number of grid points
L = 5000;        % Length of the spatial domain
dx = L / (Nx-1); % Spatial step size
x = linspace(0, L, Nx); % Spatial points
A=0.1963;
f = 0.008;
D = 500e-3;
rho=1.2;
c = 348.5;
P_atm = 1.013e5; % Atmospheric pressure in Pa
tspan = [0, 3600]; % Time span
tcut = [tspan(1) 600 1800 tspan(end)];
mcut = [0 509.28 0];

options = odeset('RelTol', 1e-6, 'AbsTol', 1e-6);

t_solution = [];
P_solution = [];
m_solution = [];

for i = 1:numel(mcut)
  tspan = [tcut(i) tcut(i+1)] ; 
  
  if i == 1
     % Initial conditions
     P0 = 5e6*ones(Nx,1); % Initial condition for pressure (P)
     m0 = zeros(Nx, 1); % Initial condition for mass flow rate (m)
     m0(end) = mcut(1);
  else
     P0 = u(end, 1:Nx).';
     m0 = u(end, Nx + 1:end).';
     m0(end) = mcut(i); 
  end
  % Initial conditions vector
  u0 = [P0;m0];
  % Solve the system of ODEs using ode15s
  [t, u] = ode15s(@(t, u) MyODEs(t, u, Nx, dx, f, c, D), tspan, u0,options);

  % Extract the solutions for P and m
  t_solution = [t_solution;t];
  P_solution = [P_solution;u(:, 1:Nx)];
  m_solution = [m_solution;u(:, Nx + 1:end)];
  
  if t(end) < tspan(2) 
      break
  end
end

%% when m =kg/s*m^2
Q_n= (m_solution./rho); 
Q=(Q_n .*A);

nodes_to_plot = 1:7:Nx;

% Plotting P_solution for selected nodes with transparency
figure;
hold on;
colors = rand(length(nodes_to_plot), 3);
for idx = 1:length(nodes_to_plot)
    i = nodes_to_plot(idx);
    plot(t_solution , P_solution(:, i), 'LineWidth', 0.3, 'DisplayName', ['Node ' num2str(i)], 'Color', colors(idx,:), 'LineStyle', '-');
end
hold off;
xlabel('Time (t)');
ylabel('Pressure (P)');
title('Pressure (P) at Selected Nodes over Time');
 legend('show');

% % Choose Node 1 (x = 1) for plotting
 node1_idx = 25;

% % Extract the pressure and Q values at Node 1
Q_node25 = Q(:, node1_idx);
figure;
plot(t_solution,  Q_node25 , 'LineWidth',2);
xlabel('Time (t)');
ylabel('Q at Node 25');
title('Q at Node 25 over Time with covective term');

function dudt = MyODEs(t, u, Nx, dx, f, c, D)
    P = u(1:Nx);
    m = u(Nx + 1:end);
    dPdt = zeros(Nx, 1);
    dmdt = zeros(Nx, 1);
    % Boundary Condition
    P(1) = 5e6; % P(1) boundary condition
    
    % Nozzle at i=25 (Algebraic relationships)
    p_hst = 503.51e3; % p_hst = rho * g * d (where d = 40m)
    Cd = 0.55;
    gamma = 1.4;
    rho1 = 1.2;
   
    % Mass flowrate of nozzle 
    m_nz =  Cd *sqrt((2 * gamma / (gamma - 1)) * P(25) * rho1 * (1 - (p_hst / P(25))^((gamma - 1) / gamma)) * (p_hst / P(25))^(2 / gamma));
    %fprintf(' m_nz: %f kg/s*m^2\n',  m_nz);

     
    % Equation
     dmdt(25) =-(2*m(25)*c^2/P(25)*((- 3*m(25)+4*m(26)-m(27)) / (2 * dx)))+(c^2*m(25)^2/P(25)^2)*((- 3*P(25)+ 4*P(26)-P(27)) / (2*dx)) - ((-3*P(25)+4*P(26)-P(27)) / (2*dx)) - (f * m(25) * abs(m(25)) * c^2) / (2 *D* (P(25))); % Forward discretization for dm/dt at i=25
     boundary_cond_m(25)=m(25)+m_nz;
     dPdt(25) =-c^2 * (3*boundary_cond_m(25) - 4*m(24) + m(23)) / (2*dx); % Backward discretization of dP/dt at i=25
    
     % Equation
     dmdt(1) =-(2*m(1)*c^2/P(1)*((- 3*m(1)+4*m(2)-m(3)) / (2 * dx)))+(c^2*m(1)^2/P(1)^2)*((- 3*P(1)+ 4*P(2)-P(3)) / (2*dx)) - ((-3*P(1)+4*P(2)-P(3)) / (2*dx)) - (f * m(1) * abs(m(1)) * c^2) / (2 *D* (P(1))); % Forward discretization for dm/dt at i=1
     dPdt(Nx) =-c^2 * (3*m(Nx) - 4*m(Nx - 1) + m(Nx-2)) / (2*dx); % Backward discretization of dP/dt at i=Nx
    
 % Central discretization for in-between grid points (i=2 to 49)
for i = 2:49
    if i == 25
        % Skip the calculation for dPdt(25) since it's already calculated outside the loop
         %dPdt(25) = ...; 
         %dmdt(25) =...;
    else
        % Calculate dP/dt and dm/dt for in-between grid points using central discretization
        dPdt(i) = -c^2 * (m(i+1) - m(i-1)) / (2 * dx); % Central discretization of dP/dt
        dmdt(i) =-2*m(i)*c^2/P(i)*((m(i + 1) - m(i - 1)) / (2 * dx))+(c^2*m(i)^2/P(i)^2)*(P(i+1) - P(i-1)) / (2*dx) - (P(i+1) - P(i-1)) / (2*dx) - (f * m(i) * abs(m(i))* c^2 ) /( 2*D*P(i)); % Central discretization of dm/dt
    end
end
    % Combine dPdt and dmdt into a single vector
    dudt = [dPdt; dmdt];
end


%% why the graph is bulge at Node 25, is it correct or wrong?

Answer 1

[Torsten]

Choose less output times from the solver, e.g. by setting

  tspan = linspace(tcut(i), tcut(i+1),(tcut(i+1)-tcut(i))/10)  

instead of

  tspan = [tcut(i) tcut(i+1)] ; 

[MAYUR MARUTI DHIRDE]

I have a question regarding numerical stability which can be seen in the graphs, At node 25 I am using backward and forward discretization. For in between point from 2 to 49 I am using central discretization. These varying discretization is causing a stability issue. Is there any solution to solve this problem.

[Torsten]

I told you my opinion about your discretization. I won't try to heal something that I think is not curable.