This code should work, but you must specify your initial conditions. You might want to look into other ode solvers, as with these specs the output looks pretty stiff.
%%%%%%%% Specify your volume and initial conditions!%%%%%%%%%%%
Vspan = [0 200];
Initial_conditions = [2, 1, 0]; %In the order of Fa0, Fb0 & Fc0
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[t, y] = ode45(@your_ode_solver_fcn,Vspan,Initial_conditions);
plot(t,y(:,1),t,y(:,2),t,y(:,3))
function out = your_ode_solver_fcn(t,in)
%Specify the input:
Fa = in(1);
Fb = in(2);
Fc = in(3);
%Define your relationships
Kc = 0.01;
Ft = Fa+ Fb+ Fc;
Co = 1;
K = 10;
Kb = 40;
ra = - (K*Co/Ft)* (Fa-Co^2*Fb*Fc^2/ (Kc*Ft^2));
Ka = 1;
Ra = Ka*Co*Fa/Ft;
Rb = Kb*Co*Fb/Ft;
%Specify the output
dFa = ra-Ra;
dFb = -ra-Rb;
dFc = -2*ra;
out = [dFa; dFb; dFc];
end
