Error using Dsolve: "Invalid indexing or function definition. Indexing must follow MATLAB indexing. Function arguments must be symbolic variables, and..."

Hello,

Upon executing the dsolve function using the code provided below, I receive the error message, "Invalid indexing or function definition. Indexing must follow MATLAB indexing. Function arguments must be symbolic variables, and function body must be sym expression." My indexing appears to be correct and I do not believe that my inputs go againt that of the function defintion, so I am unsure why this error message is given.

syms L(t) S(t) B(t)
total=10000;
f1=20;
f2=10;
f3=0.00095;
r1=0.1959;
r2=0.095;
r3=0.110;
c1=rand;
c2=rand;
c3=rand;
M1=rand;
M2=rand;
M3=rand;
gamma1=0.0001;
gamma2=0.0001;
gamma3=0.01;
k1=0.02;
k2=0.02;
k3=0;

    ode1 = diff(L) == f1/(f1+f2+f3)*k1*L + f1/(f1+f2+f3)*k2*S + r1*L*(1-(L/M1)) - k1*L - (c1*L/(1+gamma1*L));
    ode2 = diff(S) == f2/(f1+f2+f3)*k1*L + f2/(f1+f2+f3)*k2*S + r2*S*(1-(S/M2)) - k2*S - (c2*S/(1+gamma2*S));
    ode3 = diff(B) == f3/(f1+f2+f3)*k1*L + f3/(f1+f2+f3)*k2*S + r3*B*(1-(B/M3)) - k3*B - (c3*B/(1+gamma3*B));
    odes = [ode1; ode2; ode3]

odes(t) = 

    cond1 = L(0) == (20/30.00095)*total;
    cond2 = S(0) == (10/30.00095)*total;
    cond3 = B(0) == (0.00095/30.00095)*total;
    conds = [cond1; cond2; cond3];

    [LSol(t),SSol(t),BSol(t)] = dsolve(odes,conds)

Warning: Unable to find symbolic solution.

Error using sym/subsindex
Invalid indexing or function definition. Indexing must follow MATLAB indexing. Function arguments must be symbolic variables, and function body must be sym expression.

Answer 1

[Star Strider]

The differential equations are nonlinear, so it is highly unlikely that there is an analytic solution for them.

Integrate them numerically instead (using the appropriate initial conditions and time span) —

syms L(t) S(t) B(t) T Y
total=10000;
f1=20;
f2=10;
f3=0.00095;
r1=0.1959;
r2=0.095;
r3=0.110;
c1=rand;
c2=rand;
c3=rand;
M1=rand;
M2=rand;
M3=rand;
gamma1=0.0001;
gamma2=0.0001;
gamma3=0.01;
k1=0.02;
k2=0.02;
k3=0;

    ode1 = diff(L) == f1/(f1+f2+f3)*k1*L + f1/(f1+f2+f3)*k2*S + r1*L*(1-(L/M1)) - k1*L - (c1*L/(1+gamma1*L));
    ode2 = diff(S) == f2/(f1+f2+f3)*k1*L + f2/(f1+f2+f3)*k2*S + r2*S*(1-(S/M2)) - k2*S - (c2*S/(1+gamma2*S));
    ode3 = diff(B) == f3/(f1+f2+f3)*k1*L + f3/(f1+f2+f3)*k2*S + r3*B*(1-(B/M3)) - k3*B - (c3*B/(1+gamma3*B));
    odes = [ode1; ode2; ode3]

odes(t) = 

    cond1 = L(0) == (20/30.00095)*total;
    cond2 = S(0) == (10/30.00095)*total;
    cond3 = B(0) == (0.00095/30.00095)*total;
    conds = [cond1; cond2; cond3];

%     [LSol(t),SSol(t),BSol(t)] = dsolve(odes,conds)

[VF,Sbs] = odeToVectorField(odes)

VF = 

Sbs = 

odes_fcn = matlabFunction(VF, 'Vars',{T,Y})

odes_fcn = function_handle with value:
    @(T,Y)[(Y(1).*3.526702293526988-1.0).*Y(1).*(-1.9e+1./2.0e+2)-Y(1).*1.333354443775947e-2+Y(2).*6.666455562240528e-3-(Y(1).*8.300657073938436e-1)./(Y(1)./1.0e+4+1.0);Y(2).*(Y(2).*1.217628709998036e+1-1.0).*(-1.959e-1)+Y(1).*1.333291112448106e-2-Y(2).*6.667088875518943e-3-(Y(2).*2.764367447602704e-1)./(Y(2)./1.0e+4+1.0);Y(1).*6.333132784128502e-7+Y(2).*6.333132784128502e-7-Y(3).*(Y(3).*1.237874710934431-1.0).*(1.1e+1./1.0e+2)-(Y(3).*6.96175398487631e-1)./(Y(3)./1.0e+2+1.0)]

[t,y] = ode15s(odes_fcn, [0 50], zeros(1,3)+0.1);

figure
plot(t,y)
grid
legend(string(Sbs), 'Location','best')

.