Modelling Assignment Function

function [solns] = ModellingAssignment1BFuns(t, T1, T2, C1, C2)

% time (s) is independent variable
%Definition of global variables
global t
global M1 % Mass in tank 1 kg
global M2 % Mass in tank 2 kg
global cp % Specific heat capacity of reactant J.kg-1 °C-1
global T1 % Temperature in tank 1 °C
global T2 % Temperature in tank 2 °C
global THi % Temperature of jacket fluid in °C
global THo % Temperature of jacket fluid out °C
global THm % Temperature of jacket between tank 1 and 2 °C
global E % Activation energy J.mol-1
global R % Ideal gas constant J.mol-1K-1
global k % Rate constant at (mol.kg-1)1-n s-1
global T1i % Initial temperature in tank 1 °C
global T2i % Initial temperature in tank 2 °C
global F % Flowrate between tanks kg.s-1
global q % Flowrate of heating fluid in jackets kg.s-1
global cpw % Specific heat capacity of heating fluid J.kg-1 °C-1
global C1 % Concentration in tank 1 mol.kg-1
global C2 % Concentration in tank 2 mol.kg-1
global Cf % Feed concentration to tank 1 mol.kg-1
global UA % Overall htc jacket to tanks W.m-2°C-1
global C1i % Initial concentration in tank 1 mol.kg-1
global C2i % Initial concentration in tank 2 mol.kg-1

% Make sure the time arrays used in the ode solver is the same interval and
% length as in the experimental data

% Define an array to hold the solutions
[solns] = zeros(1, 5); % 1 row, 4 cols (time, temp1/temp2, conc1/conc2)
solns(1) = t;
% These are the ODEs that model the system
% 𝑀1𝑐𝑝*𝑑𝜃/𝑑𝑡 = 𝐹𝑐_𝑝(𝜃𝐹 − 𝜃1) + 𝑈1𝐴1∆𝜃𝑙𝑚1 𝑓𝑜𝑟 𝑡 > 0
% 𝑀2𝑐𝑝*𝑑𝜃/𝑑𝑡 = 𝐹𝑐𝑝(𝜃1 − 𝜃2) + 𝑈2𝐴2∆𝜃𝑙𝑚2 𝑓𝑜𝑟 𝑡 > 0
% 𝑀1𝑑𝐶1/𝑑𝑡 = 𝐹(𝐶𝐹 − 𝐶1) − 𝑘1𝐶𝑀1 𝑓𝑜𝑟 𝑡 > 0
% 𝑀2𝑑𝐶2/𝑑𝑡 = 𝐹(𝐶1 − 𝐶2) − 𝑘2𝐶𝑀2 𝑓𝑜𝑟 𝑡 > 0

% Algebraic equations
% Arrhenius law
ko = 10;
k = ko*exp(-E/(R *(T1 + 273.15)));
% Calculate THo and THm
THo = T1 + (THm - T1)*exp(-UA/F/cpw);
THm = T2 + (THi - T2)*exp(-UA/F/cpw);

% Calculate log mean temp differences
lm1 = (THm - THo)/log((THm - T1)/(THo - T1));
lm2 = (THi - THm)/log((THi - T2)/(THm - T2));

% ODEs
% Calculate temperature for tank 1
dT1_dt = (F*cp*(T - T1i) + UA*lm1)/ (M1*cpw);
solns(2) = dT1_dt;
% Calculate temperature for tank 2
dT2_dt = (F*cp*(Tf - T2i) + UA*lm2)/ (M2*cpw);
solns(3) = dT2_dt;

% Calculate concentration in tank 1
dC1_dt = (F(Cf - C1) - k*C1*M1) / M1;
solns(4) = dC1_dt;
% Calculate concentration in tank 2
dC2_dt = (F(C1 - C2) - k*C2*M2) / M2;
solns(5) = dC2_dt;