一个论文复现，不太会matlab，有大佬帮忙编一下代码么，用ai写的代码跑不出来

%% TNT密闭燃烧-不同泄压直径压力对比（按照论文公式计算）

%% 1. 核心物理参数（按表3.1设置）

rho_e = 1620;           % 炸药密度 [kg/m³]

R = 337;                % 气体常数 [J/(kg·K)]

chi = 0.98;             % 热量修正系数

varphi1 = 0.95;         % 流量修正系数

gamma = 1.3;            % 比热比

D_id = 0.0025;            % 炸药厚度 [m]

R_ex = 0.025;           % 光斑孔直径 [m]（25mm直径）


% TNT燃烧特性（论文参数）

a = 0.009 * 1e-2;  % 转换为 m/s (原单位 cm/s)

b = 0.05 * 1e-2;   % 转换为 m/(Pa·s) (原单位 cm/(Pa·s))

n = 1;

% 压力指数


% 固定燃烧温度（按表3.1）

T = 3400;               % 燃料燃烧温度 [K]，固定值


% 初始条件（论文设定）

p0 = 2e6;               % 初始压力 [Pa] → 2MPa

rho0 = p0 / (R * T);    % 通过状态方程计算初始气体密度 [kg/m³]

e0 = 0;                 % 初始已燃厚度 [m]

V0 = pi * (R_ex/2)^2 * D_id; % 初始燃烧腔体积 [m³]


% 初始状态向量 [V1; e; rho1; p1]

y0 = [V0; e0; rho0; p0];


% 定义不同泄压孔直径（单位：m），可根据需要增减

diameters = [0.002, 0.005, 0.010, 0.08]; % 2mm、5mm、10mm、15mm

colors = ['r', 'g', 'b', 'm'];            % 对应曲线颜色

labels = {'2mm', '5mm', '10mm', '15mm'};  % 对应图例标签


%% 2. 循环计算不同泄压直径下的压力

figure('Color','w','Position',[100,100,900,500]);

hold on; grid on;


% 存储结果用于后续分析

results = cell(length(diameters), 1);


for i = 1:length(diameters)

% 计算当前直径对应的排气面积

    A_t = pi * (diameters(i)/2)^2;


% 求解微分方程

    tspan = [0, 0.1]; % 延长计算时间到10秒

    options = odeset('RelTol',1e-6,'AbsTol',1e-8,'NonNegative',[1,2,3,4]);


try

        [t, y] = ode45(@(t,y) combustion_ode_paper(t,y,rho_e,R,chi,varphi1,gamma,D_id,R_ex,A_t,a,b,n,T), ...

            tspan, y0, options);


% 提取压力并绘图（转换为MPa）

        p = y(:,4)/1e6;

        plot(t, p, ['-', colors(i)], 'LineWidth',1.5, 'DisplayName', labels{i});


% 保存结果

        results{i} = struct('t', t, 'p', p, 'diameter', diameters(i));


catch ME

        warning('直径 %.3f m 计算失败: %s', diameters(i), ME.message);

end

end


%% 3. 绘图美化

xlabel('时间 (s)','FontSize',12);

ylabel('压力 (MPa)','FontSize',12);

title('不同泄压孔直径下TNT密闭燃烧压力变化曲线（论文公式）','FontSize',14);

legend('Location','best','FontSize',10);

ylim([0, 100]);  % 对应论文90MPa安全阈值

xlim([0, 10]);   % 限制x轴范围


% 添加安全阈值线

yline(90, 'k--', 'LineWidth', 1.5, 'DisplayName', '安全阈值(90MPa)');


set(gca,'FontSize',10);

grid on;


%% 4. 分析结果：判断是否超过安全阈值

fprintf('\n=== 安全阈值分析 (90MPa) ===\n');

for i = 1:length(results)

if ~isempty(results{i})

        max_p = max(results{i}.p);

        is_safe = max_p < 90;

        status = '安全';

if ~is_safe

            status = '危险(可能爆炸)';

end

        fprintf('直径 %.1f mm: 最大压力 = %.2f MPa, 状态: %s\n', ...

            diameters(i)*1000, max_p, status);

end

end


%% 6. 可选：绘制燃烧面积和已燃厚度随时间变化

if exist('results', 'var') && ~isempty(results)

    figure('Color','w','Position',[100,100,1200,400]);


% 子图1：燃烧面积变化

    subplot(1,2,1);

    hold on; grid on;

for i = 1:length(results)

if ~isempty(results{i})

% 重新计算燃烧面积（为了绘图）

            t = results{i}.t;

            e = interp1(results{1}.t, results{1}.y(:,2), t); % 近似使用第一个结果的厚度

            A_b = pi * D_id * (R_ex + 2*e);

            plot(t, A_b, ['-', colors(i)], 'LineWidth',1.5, 'DisplayName', labels{i});

end

end

    xlabel('时间 (s)','FontSize',12);

    ylabel('燃烧面积 (m²)','FontSize',12);

    title('燃烧面积随时间变化','FontSize',14);

    legend('Location','best','FontSize',10);


% 子图2：已燃厚度变化

    subplot(1,2,2);

    hold on; grid on;

for i = 1:length(results)

if ~isempty(results{i})

            t = results{i}.t;

            e = interp1(results{1}.t, results{1}.y(:,2), t); % 近似使用第一个结果的厚度

            plot(t, e*1000, ['-', colors(i)], 'LineWidth',1.5, 'DisplayName', labels{i}); % 转换为mm

end

end

    xlabel('时间 (s)','FontSize',12);

    ylabel('已燃厚度 (mm)','FontSize',12);

    title('已燃厚度随时间变化','FontSize',14);

    legend('Location','best','FontSize',10);

end


%% 5. 微分方程函数定义（按照论文公式）

function dydt = combustion_ode_paper(~, y, rho_e, R, chi, varphi1, gamma, D_id, R_ex, A_t, a, b, n, T)

% 状态变量

    V1 = y(1);  % 燃烧腔体积 [m³]

    e = y(2);   % 已燃厚度 [m]

    rho1 = y(3);% 腔内气体密度 [kg/m³]

    p1 = y(4);  % 腔内压力 [Pa]


% 1. 燃烧面积（柱状燃烧，按论文公式）

% A_b = (R_ex + 2e) * π * D_id

    A_b = pi * D_id * (R_ex + 2*e);


% 2. 燃烧速度（指数燃烧定律）

    u_e = a + b * p1^n;


% 3. 质量流入（燃烧产生）

    m_dot_in = rho_e * A_b * u_e;


% 4. 质量流出（壅塞流量，按论文公式(3-24)形式）

% 注意：论文公式(3-24)中的质量流出项为：

% m_dot_out = p1 * ((1+gamma)/2) * (gamma*R*T/(gamma+1))^(1/2) * φ1 * A_t

% 但公式(3-27)中使用了不同的形式，这里采用更完整的形式


% 计算壅塞流量系数

% 根据公式(3-27)中的质量流出相关项

    term1 = sqrt(gamma / (R * T));

    term2 = (2 / (gamma + 1))^((gamma + 1) / (2 * (gamma - 1)));

    m_dot_out = varphi1 * A_t * p1 * term1 * term2; % 添加流量修正系数


% 5. 质量守恒方程（论文公式3-24）

% V1 * dρ1/dt = ρ_e * A_b * u_e - m_dot_out

    drho1_dt = (rho_e * A_b * u_e - m_dot_out) / V1;


% 6. 能量守恒方程（论文公式3-27）

% V1 * dp1/dt = ρ_e * A_b * u_e * R * χ * T * γ - 

%               p1 * ((γ+1)/2) * (1/(γ-1)) * sqrt(2γRT/(γ+1)) * φ1 * γ * R * T * A_t

% 能量流入项

    energy_in =  m_dot_in * R * chi * T * gamma;


% 能量流出项（按公式3-27）

    energy_out_term = p1 * ((gamma+1)/2) * (1/(gamma-1)) * ...

                     sqrt(2*gamma*R*T/(gamma+1)) * varphi1 * gamma * R * T * A_t;


% 压力变化率

    dp1_dt = (energy_in - energy_out_term) / V1;


% 7. 体积变化率（论文公式3-19）

    dV1_dt = A_b * u_e;


% 8. 已燃厚度变化率（论文公式3-20）

    de_dt = u_e;


% 返回导数向量

    dydt = [dV1_dt; de_dt; drho1_dt; dp1_dt];


% 物理量保护：防止出现负值或不合理的值

% 确保密度和压力非负

    dydt(3) = max(dydt(3), -rho1/1e-3); % 限制密度变化率

    dydt(4) = max(dydt(4), -p1/1e-3);   % 限制压力变化率


% 确保体积和已燃厚度非负且增长合理

    dydt(1) = max(dydt(1), 0);

    dydt(2) = max(dydt(2), 0);

end