I am trying to write the MATLAB code of the question ofresearch paper Drive link of paper-http​s://drive.​google.com​/file/d/12​N5ugCymHKB​ubfqOlM2LU​3DVBxdrCoU​R/view?usp​=drive_lin​k

function batch_distillation()
% Binary Data
ki = [1.51; 1.13; 0.95; 0.91];
alphaij = [5.91, 2.51, 1.29, 1.00;
    1/5.91, 1, 1/2.51, 1/1.29;
    1/2.51, 1/2.51, 1, 1/1.29;
    1/1.29, 1/1.29, 1/1.29, 1];

% Given data
F_feed = 400; % kmol/day
F_i = [100; 100; 100; 100]; % kmol of each component in the feed
x_i = [0.97; 0.97; 0.97; 0.97]; % molar fractions in the feed
K = ki .* x_i; % Equilibrium constants
alpha = alphaij .* (x_i ./ x_i'); % Relative volatilities

% Time parameters
t_span = [0 10]; % Set the time span for batch distillation (adjust as needed)
t_steps = 100; % Number of time steps

% Initial conditions
L0 = 0; % Initial liquid holdup
V0 = 0; % Initial vapor holdup
z_L0 = x_i; % Initial liquid composition (same as feed composition)
z_V0 = x_i; % Initial vapor composition (same as feed composition)

% Solve the batch distillation using ODE solver
options = odeset('RelTol', 1e-6, 'AbsTol', 1e-9);
[t, y] = ode45(@(t,y) batch_distillation_ode(t, y, F_feed, ki, alphaij, x_i, t_span), t_span, [L0; V0; z_L0; z_V0], options);

% Extract results
L = y(:, 1);
V = y(:, 2);
z_L = y(:, 3:6);
z_V = y(:, 7:10);

% Plot results
figure;
plot(t, L, 'b', 'LineWidth', 2);
hold on;
plot(t, V, 'r', 'LineWidth', 2);
xlabel('Time (hours)');
ylabel('Moles (kmol)');
legend('Liquid holdup (L)', 'Vapor holdup (V)');
title('Moles vs Time');

figure;
plot(t, z_L, 'LineWidth', 2);
hold on;
plot(t, z_V, 'LineWidth', 2);
xlabel('Time (hours)');
ylabel('Molar Fraction');
legend('Benzene (L)', 'Toluene (L)', 'Ethylbenzene (L)', 'o-xylene (L)',...
    'Benzene (V)', 'Toluene (V)', 'Ethylbenzene (V)', 'o-xylene (V)');
title('Molar Fractions vs Time');
end

function dydt = batch_distillation_ode(t, y, F_feed, ki, alphaij, x_i, t_span)
% Unpack variables
L = y(1);
V = y(2);
z_L = y(3:6);
z_V = y(7:10);

% Calculate liquid and vapor flow rates
F_L = L / t_span(end);
F_V = V / t_span(end);

% Calculate liquid and vapor compositions
x_L = z_L ./ sum(z_L);
x_V = z_V ./ sum(z_V);

% Calculate liquid and vapor phase rates
r_L = F_L .* x_L;
r_V = F_V .* x_V;

% Calculate material balances
dLdt = F_feed - F_L;
dVdt = F_L - F_V;
dz_Ldt = (r_L - z_L) ./ L;
dz_Vdt = (r_V - z_V) ./ V;

% Pack the derivatives
dydt = [dLdt; dVdt; dz_Ldt; dz_Vdt];
end

[nayan gupta]

please help in correcting this code.

it is some error in some lines

[Walter Roberson]

what error message are you getting?

[Aditya]

It seems that due to your initialization all the values are coming to be INF for z_L and z_V, please do check your initialization. if the issue still persists please do share the paper that you are refering to.