I need matlab code to do this, I have tried 3 times already with my code that doesnt completely work underneath

Normal butane is to be isomerised to isobutane in a PFR. For this example, the reaction is elementary and not reversible (ie neglect equilibrium). The new specific reaction rate constant of 10 h-1 was re-measured at 380K after a students data was found to be erroneous. The feed for this reactor will enter at 330K and will use an existing reactor of 6m3. It will process 200kmol/hr of a mixture of 40% butane and 60% N2.

A) Plot the conversion X and reactor Temperature along the reactor assuming adiabatic conditions and plot the conversion X and reactor Temperature along the reactor using an assumed constant ambient temperature Ta of 25C and a overall heat transfer constant Ua of 2 kJ/h.K

B) Now conduct a hot spot sensitivity analysis for ambient temperatures of 10, 20, 30, 40, 50C to determine the hottest temperature in the reactor and where it occurs. Also determine the overall conversion at the end of the reactor in each case.

My code

% Constants

deltaHrxn = -9000; % J/mol

Ea = 50000; % kJ/mol

Cp_butane = 141; % J/mol.K

Cp_i_butane = 151; % J/mol.K

Cp_N2 = 29.14; % J/mol.K

deltaCp = 10; % J/mol.K

Butane0 = 15; % kmol/m^3

FeedFlowRate = 200; % kmol/hr

ButaneFraction = 0.4;

N2Fraction = 0.6;

% Convert temperatures to Kelvin

To = 330;

% Solve the ODEs

vspan = 0:0.01:6; % R-Vol varying from 0 to 6 m^3

IC = [To, 0]; % Initial conditions [Temperature, Conversion]

[vsol, varsol] = ode23(@od_sys, vspan, IC);

% Extract variables from the solution

Temperature = varsol(:, 1);

Conversion = varsol(:, 2);

% Plot conversion and reactor temperature along the reactor (adiabatic conditions)

figure;

subplot(2, 1, 1);

plot(vsol, Conversion, 'm');

title('Conversion of n-Butane along the Reactor (Adiabatic Conditions)');

xlabel('Volume (m^3)');

ylabel('Conversion X');

subplot(2, 1, 2);

plot(vsol, Temperature, 'r');

title('Reactor Temperature along the Reactor (Adiabatic Conditions)');

xlabel('Volume (m^3)');

ylabel('Temperature (K)');

% Modify ambient temperature and overall heat transfer constant for constant ambient temperature case

Ta = 25 + 273.15; % Constant ambient temperature of 25C

Ua = 2; % kJ/h.K

% Solve the ODEs with modified parameters

[vsol1, varsol] = ode23(@od_sys, vspan, IC);

% Extract variables from the solution

Temperature = varsol(:, 1);

Conversion = varsol(:, 2);

% Plot conversion and reactor temperature along the reactor (constant ambient temperature)

figure;

subplot(2, 1, 1);

plot(vsol, Conversion, 'm');

title('Conversion of n-Butane along the Reactor (Constant Ambient Temperature)');

xlabel('Volume (m^3)');

ylabel('Conversion X');

subplot(2, 1, 2);

plot(vsol, Temperature, 'r');

title('Reactor Temperature along the Reactor (Constant Ambient Temperature)');

xlabel('Volume (m^3)');

ylabel('Temperature (K)');

% Sensitivity analysis for ambient temperatures

ambientTemps = [10, 20, 30, 40, 50]; % Celsius

hottestTemps = zeros(size(ambientTemps));

overallConversions = zeros(size(ambientTemps));

function contdiff = od_sys(T,var)

T=var(1); % Input of temperature (K)

X=var(2); % Input of Conversion

Ua=2*10^3;

Fao= 0.4*200000;% Flow rate of Component mol/hr

Cao=15000;% Initial Concentration of A (mol/m^3)

dCp=151-141;

k= dCp*exp((50000/8.314)*((1/380)-(1/T))); % Rate constant at any temperature T

ra= -k*Cao*(1-X); % Rate of Reaction of Production of A i.e. why it is -ve

% Because A is consumed in the reaction.

delH=-9000; % Heat of Reaction at any temp.( Since del Cp is 0)

% dT/dV (Energy Balnce)

contdiff(1,1)= ((Ua*(273.15+50-T))+(ra*delH))/(Fao*(184.71+dCp*X));

%dX/dV (Mole Balance)

contdiff(2,1)=-ra/Fao;

end

%__________ Script File_____________