Error using ode45

Hello. I'm trying to reproduce this model: http://www.ual.es/~jfernand/ProcMicro70801207/tema-3---cinetica-del-crecimiento-de-microalgas/3-7-optimizacion-batch.html This is my code:

clear all; clc;
Io=0.00150;
k=0.075;
mumax=0.045;
kappa=0.075;
a=0.00018;
Leq=0.05;
ka=0.19;
Cb=100;
m=0.005;

Iav = (Io/(ka*Cb*Leq))*(1-exp(-ka*Cb*Leq));

mo= mumax * (Iav/(2*a)) * ( 1+k+(a/Iav)- sqrt(( 1-k-(a/Iav))^2 + (4*k) ));

f=@(t,Cb) Cb(mo-m);

tf=1000;
tspan=[0 tf];
Cb0=50;

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

plot(t,Cb);

But I get the next error:

Attempted to access Cb(0.0362687); index must be a positive integer or logical.

Error in @(t,Cb)Cb(mo-m)

Error in odearguments (line 88)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 114)
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in OptFBPlacas (line 22)
[t,Cb]=ode45(f,tspan,Cb0);

What am I doing wrong? Help me, please.

Answer 1

[Star Strider]

You’re doing implied multiplication (that MATLAB does not recognise), so it thinks you are addressing an array, and is throwing the error. I cannot read the reference you posted, so I also have to

See if this works:

Iav = @(Cb) (Io./(ka.*Cb.*Leq)).*(1-exp(-ka.*Cb.*Leq));

mo = @(Cb) mumax .* (Iav(Cb)./(2*a)) .* ( 1+k+(a./Iav(Cb))- sqrt(( 1-k-(a./Iav(Cb))).^2 + (4*k) ));

f=@(t,Cb) Cb.*(mo(Cb)-m);

Note that in the MathCad example, both ‘Iav’ and ‘mu’ are functions of several variables.

Your constants span a few orders of magnitude, so if ode45 gets stuck or gives strange results, see if ode15s will work.

Note — This is UNTESTED CODE so you may have to experiment with it to get it to work. My changes should at least solve some of your problems, and the revised code will not throw the error it did earlier.

[Mario Gonzaléz Sandoval]

Thank you! :D That was the problem.

[Star Strider]

My pleasure!

I aplolgise for the delay — my laptop crashed.