Untitled

function run_simulation
    clear;
    close all;
    clc;

    global F V C0 T0 Cp Rho k0 E R dH Th Cph Rhoh Fh
    F = 1;
    V = 5;
    C0 = 5;
    T0 = 300;
    Cp = 1;
    Rho = 1e3;
    k0 = 1e7;
    E = 1.104e4;
    R = 1.987;
    dH = 1.1e5;
    Th = 363;
    Cph = 1;
    Rhoh = 1e3;
    Fh = 20;

    t_start = linspace(0, 20, 1000);
    t_end = linspace(20, 40, 1000);
    init_vals = [C0; T0];

    [t1, y1] = ode15s(@model_eqns, t_start, init_vals);
    
    Fh = 40;
    init_vals_2 = y1(end,:);
    [t2, y2] = ode15s(@model_eqns, t_end, init_vals_2);

    t_combined = [t1; t2];
    y_combined = [y1; y2];

    figure('Position', [100 100 800 600]);

    subplot(2,1,1);
    plot(t_combined, y_combined(:,1), 'LineWidth', 2);
    grid on;
    xlabel('Time (min)', 'FontSize', 12);
    ylabel('Concentration (kmol/m^3)', 'FontSize', 12);
    title('Concentration Profile', 'FontSize', 14);
    hold on;
    plot([20 20], [2 5.5], '--k', 'LineWidth', 1.5);
    text(21, 4, 'Flow rate change', 'FontSize', 10);
    xlim([0 40]);
    ylim([2 5.5]);

    subplot(2,1,2);
    plot(t_combined, y_combined(:,2), 'LineWidth', 2);
    grid on;
    xlabel('Time (min)', 'FontSize', 12);
    ylabel('Temperature (K)', 'FontSize', 12);
    title('Temperature Profile', 'FontSize', 14);
    hold on;
    plot([20 20], [295 335], '--k', 'LineWidth', 1.5);
    text(21, 315, 'Flow rate change', 'FontSize', 10);
    xlim([0 40]);
    ylim([295 335]);
end

function dydt = model_eqns(t, y)
    global F V C0 T0 Cp Rho k0 E R dH Th Cph Rhoh Fh
    
    CA = y(1);
    T = y(2);
    
    k = k0 * exp(-E/(R*T));
    r = k * CA;
    
    Q = 12 * Fh^2 * (Th - T);
    
    dCAdt = (F/V) * (C0 - CA) - r;
    dTdt = (F/V) * (T0 - T) + Q / (Rho * V * Cp) + (-dH) * r / (Rho * Cp);
    
    dydt = [dCAdt; dTdt];
end