Solving a 4 ODE system in MATLAB using ode45

I am not very used to MATLAB and I'm trying to solve the following problem using MATLAB ode45, however, it's not working.

I was working on a problem in reaction engineering, using a Semi-Batch Reactor. The reaction is given by A + B ---> C + D A is placed in the reactor and B is being continuously added into the reactor with a flowrate of v0 = 0.05 L/s. Initial volume is V0 = 5 L. The reaction is elementary. The reaction constant is k = 2.2 L/mol.s. Initial Concentrations: for A: 0.05 M, for B: 0.025 M.

Performing a mole balance of each species in the reactor, I got the following 4 ODEs, and the expression of V (volume of the reactor is constantly increasing)

Solving this system and plotting the solution against time, I should get this

Note that plots of C(C) and C(D) are the same. And let's set tau = v0/V.

Now for the MATLAB code part.

I have searched extensively online, and from what I've learned, I came up with the following code.

First, I wrote the code for the ODE system

function f = ODEsystem(t, y, tau, ra, y0)
f = zeros(4, 1);
f(1) = ra - tau*y(1);
f(2) = ra + tau*(y0(2) - y(2));
f(3) = -ra - tau*y(3);
f(4) = -ra - tau*y(4);
end

Then, in the command window,

t = [0:0.01:5];
v0 = 0.05;
V0 = 5;
k = 2.2;
V = V0 + v0*t;
tau = v0./V;
syms y(t);
ra = -k*y(1)*y(2);
y0 = [0.05 0.025 0 0];
[t, y] = ode45(@ODEsystem(t, y, tau, ra, y0), t, y0); 
plot(t, y); 

However, I get this...

Please if anyone could help me fix my code. This is really annoying :)

ra should not be passed as parameter but be computed inside the ODE system. V is likewise not a constant. Symbolic expressions should be used for formula transformations, not for numerical methods. One would also have to explicitly evaluate the symbolic expression at the wanted numerical values.

function f = ODEsystem(t, y, k, v0, V0, cB0)
f = zeros(4, 1);
ra = -k*y(1)*y(2);
tau = v0/(V0+t*v0);
f(1) = ra - tau*y(1);
f(2) = ra + tau*(cB0 - y(2));
f(3) = -ra - tau*y(3);
f(4) = -ra - tau*y(4);
end

Then use the time span of the graphic, start with all concentrations zero except for A, use the concentration B only for the inflow.

t = [0:1:500];
v0 = 0.05;
V0 = 5;
k = 2.2;
cB0 = 0.025;
y0 = [0.05 0 0 0];
[t, y] = ode45(@(t,y) ODEsystem(t, y, k, v0, V0, cB0), t, y0); 
plot(t, y); 

and get a good reproduction of the reference image