Variable inside an ode

I am trying to solve this differential equation: dCb/dt = ((A*D)/(DV*h))*(Cs-Cb). I specify in Matlab the function:

f = @(t,Cb) ((A*D)/(h*DV))*(Cs-Cb)

I use ode45 to solve it:

[t,Cb] = ode45(f, tspan, Cb0);

Full code is below.

However, A is a variable and varies with time itself. It is actually A(t) =0.00155+(0.06043-0.0015)/(1+(t/4481.23119).^0.85949)

How do I add this A function to the below code?

Thanks

tspan = [0 9e4];
dv50=189/1000000000; %m, 
dv90=795/1000000000;
r0=dv90/2; %m dv90
rho=1370; %kg/m3
W=18/1000000 %kg
N= W/((4/3)*pi*r0^3*rho);
A=(4*pi*r0^2*N); %m2
MW=402.88
D= 9.9e-9*MW^-0.453  %m2/s % from D(m2/s)=9.9e-9*MW^-0.453
h= r0 %m - 
Cs= 0.032 %kg/m3
DV=0.0009 % m3, volume
Cb0=0;
f = @(t,Cb) A*D/(h*DV)*(Cs-Cb)
[t,Cb] = ode45(f, tspan, Cb0);
plot(t,Cb);
title('Bulk Concentration vs Time');
xlabel('Time (s)');
ylabel('Cb (kg/m3)');
legend('Nernst Brunner','Cb_exp');

[Adam]

Can you not just replace A in the definition of f, with its above representation in terms of t which is already an input to f?

[Nora Rafael]

This seemed to work with no errors. But I am concerned as when I plot just this A function alone I get a strange curve.

This is the equation:

and I implement it this way

SA=0.00155+(0.06043-0.0015)/(1+(t/4481.23119).^0.85949)

This is the curve that I get

This is a strange curve as this function should output an exponential decline.

I tested this function separately

function SA_=SurfaceArea(t)
SA_=0.00155+(0.06043-0.0015)/(1+(t/4481.23119).^0.85949);
end

>> t=1:9e4;

>> SA_=SurfaceArea(t)

and I get this error message:

Error using /

Matrix dimensions must agree.

Error in SurfaceArea (line 2)

SA_=0.00155+(0.06043-0.0015)/(1+(t/4481.23119).^0.85949);

Answer 1

[Star Strider]

Try this:

SA = @(t) 0.00155+(0.06043-0.0015)./(1+(t/4481.23119).^0.85949);
f = @(t,Cb) SA(t)*D/(h*DV)*(Cs-Cb)

The rest of your code is unchanged.

[Nora Rafael]

Thanks.

This seems to also work and outputted the same result as what Adam suggested. But when I come to plot SA with time, I get an error, see the last line:

[Star Strider]

My pleasure.

Note that ‘SA’ is a function, so:

figure
plot(t, SA(t))

will do what you want.

[Nora Rafael]

Super!

[Star Strider]

Thank you!