The way you've written the equations above only gives you 2 differential equations, in your odethreevariable function you code it as if the first equation was:
Either your odethreevariables function or your equations are correct, you'll have to decide which it is.
To solve this type of problems you can do something like:
K_De_0 = [1 1]; % Or a more appropriate start-guess
K_De_optimal = fminsearch(@(K_De) sum((C_B_obs-C_B_fcn(tC,K_De)).^2),K_De_0);
That is a standard least-square fit of your function C_B_fcn with parameters K_De evaluated at times tC (that should be the instances in time for which you have your estimates of C_B). The only thing you'll need to do is to integrate your ODEs inside C_B_fcn to get the C_B-curve, something like this:
function C_B = C_B_fcn(tC,K_De)
ICs=zeros(3,1);
ICs(1,1) = 0.0035;
ICs(2,1) = 0;
[t,X]=ode45(@(t,X) ode2variable(t,X,K_De),tC,ICs);
C_B = X(:,2);
end
function [dX_dt]=ode2variable(t,X,K_De)
% EDIT:
K = K_De(1);
De = K_De(2);
rc = X(1);
Cb = X(1);
rg = 0.0035;
Y0 = 182;
a = 190.47;
Jb = 1/((rg^2/(rc^2*(K*Y0)))+((rg*(rg-rc))/(rc*De)));
dX_dt = [-((Jb*(rg^2))/(Y0*rc^2)) ; Jb*a]; % ODE (system of 2 diferential equations)
end
This is not tested but should at least give you a starting point.
HTH
