Write a function
function dy = f(t,y)
that accepts t and a vector y of length 2*n with
y = [c_1;...;c_n,cbar_1;...;cbar_n]
and returns
dy = [dc_1/dt;...;dc_n/dt;dcbar_1/dt;...,dcbar_n/dt]
Then write the Runge-Kutta-method of order 4 for a system of differential equations of the form
y' = f(t,y)
and apply it to your above f.
tstart = 0.0;
tend = 5000.0;
h = (tend - tstart)/1000000;
T = (tstart:h:tend).';
n = 10; % Number of stages/bed divisions
Y0 = zeros(1,2*n);
Y0(n+1) = 9;
Y = runge_kutta_RK4(@(t,y)f(t,y,n),T,Y0);
plot(T,Y(:,10)) % Runge Kutta solution
grid on
function dy = f(t,y,n)
e = 0.38; % Bed porosity
rho = 126; % Solvent (CO2) density
W = 0.95; % CO2 flow rate, kg/h
V = 0.0004; % Extractor volume, m^3
Di = 6 * 10^-13; % Diffusion coefficient
ti = 1736.111; % Internal diffusion time
c_in = @(t) 0;
c = y(1:n);
cbar = y(n+1:2*n);
dc = zeros(1,n);
dcbar = zeros(1,n);
dcbar = -1/ti * (cbar-c/0.2);
dc(1) = -(W/rho*(c(1)-c_in(t))+(1-e)*V*dcbar(1))/(e*V);
dc(2:n) = -(W/rho*(c(2:n)-c(1:n-1))+(1-e)*V./(2:n).*dcbar(2:n))./(e*V./(2:n));
dy = [dc,dcbar];
end
function Y = runge_kutta_RK4(f,T,Y0)
N = numel(T);
n = numel(Y0);
Y = zeros(N,n);
Y(1,:) = Y0;
for i = 2:N
t = T(i-1);
y = Y(i-1,:);
h = T(i) - T(i-1);
k0 = f(t,y);
k1 = f(t+0.5*h,y+k0*0.5*h);
k2 = f(t+0.5*h,y+k1*0.5*h);
k3 = f(t+h,y+k2*h);
Y(i,:) = y + h/6*(k0+2*k1+2*k2+k3);
end
end
