Error in ode45 (line 115) and odearguments (line 90)

I'm struggling to find out what is preventing the script to not run.

These are my codes (function and run):

FUNCTION:

function dcdt = groupproject(~,c)
global umaxa
global Ksa
global umaxm
global Ksm
global qmaxa
global qmaxm
global Km
global Y
global V
global m
global F
global Mtotal
global Vol
global T
global R
global Emin
global Emax
global Rmin
global Rmax
global Kr
global Rext
global xa
global qa
global qm
global D
global So
global Kda
global Kdm

%Anode rate equations
mua = umaxa * (c(1) / (Ksa + c(1))) * (c(4) / (Km + c(4)));  %Anodophillic
qa = qmaxa *  (c(1) / (Ksa + c(1))) * (c(4) / (Km + c(4))); %Anodophillic
mum = umaxm * (c(1) / (Ksm+c(1))); %methanogenic
qm = qmaxm * (c(1) / (Ksm+c(1))); %methanogenic
%Cathode (Intracellular??) rate equations
Mred = Mtotal -1 * c(4);
IMFC = ((Emin + (Emax -1 * Emin) * exp(-1/(Kr*c(2)))) - 1 * (R*T/F)*log(Mtotal/Mred))/(Rext - 1 * Rmin + (Rmax-Rmin)*exp(-Kr*c(2)));
%Ordinary Differential Equations
dcdt = zeros(4,1);
%Anode soluble components
dcdt(1)=(-qa*xa)-(qm*.c(3))+(D*(So-.c(1)));  %Substrate
dcdt(2) = (mua * .c(2)) -1*(Kda*.c(2)) -1*(alpha_a*D*.c(2)); %anodophillic
dcdt(3) = (mum * .c(3)) -1*(Kdm*.c(3)) -1*(alpha_m*D*.c(3)); %methanogenic
%Cathode (Intracellular) Soluble components
dcdt(4) = (-1 * Y * qa) + ((V * (IMFC/(m * F))) * (1/(Vol* .c(2))));

RUN:

function rungroupproject
clear
clc
global umaxa
global Ksa
global umaxm
global Ksm
global qmaxa
global qmaxm
global Km
global Y
global V
global m
global F
global Mtotal
global YCH4
global Vol
global T
global R
global Emin
global Emax
global Rmin
global Rmax
global Kr
global Rext
global tf
global D
global Kda
global Kdm
global S

%system parameters
umaxa = 1.97; %d-1 maximum anodophillic growth rate
Ksa = 20; %mg-S L-1 Half rate constant of anodophilics
umaxm = 0.1; %d-1 maximum methanogenic growth rate
Ksm = 80; %mg-S L-1 Half rate constant of methanogens
qmaxa = 8.48; %mg-S mg-x−1 d−1 Maximum anodophilic reaction rate
qmaxm = 8.20; %mg-S mg-x−1 d−1 Maximum methanogenic reaction rate
Mtotal = 0.05; %mg-M mg-x−1 Mediator fraction
Km = 0.2 * Mtotal; %mg-M L−1 Mediator half-rate constant
Y = 22.75; %mg-M mg-S−1 Yield in Eq. (12)
V = 6634000; %mg-M molmediator−1 Mediator molar mass
m = 2; %mole−1 molmediator−1 Electrons transferred per mol of mediator
F = 96485; %A s mole−1 Faraday Constant
YCH4 = 0.3; %mL-CH4 mg-S−1 Methane Yield
Vol = 50; % mL MFC Volume
T = 298.15; %K MFC temperature
R = 8.314472; %J K−1 mol−1 Ideal gas constant
Emin = 0.1; %Volts Minimum Eocv
Emax = 0.61;
Rmin = 25;
Rmax = 2000;
Kr = 0.006;
Rext = 10;
tf=50;
Kda=0.02*umaxa;
Kdm=0.02*umaxm;
D=1;
%initial conditions
c = zeros(4,1);
c(1) = 0;  %Substrate
c(2) = 0;  %andophillics
c(3) = 0;  %methanogens
c(4) = 0;  %intracellular
K = [1 1 1 1 1 1 1 1 1 1 1 1];
Rate = [];
change = [];
%Global Options
options = odeset('MaxStep',0.1,'NonNegative',(1:4));
%Calling the Solver
S = ode45(@groupproject,(0:0.1:tf),c,options);
Rate = [Rate S.y(:,4)];
%1.5 times k
for i = 1:12
   K(i) = 1.5;
   %Kinetic Parameters
   Ksa = 20; %mg-S L-1 Half rate constant of anodophilics
   Ksm = 80; %mg-S L-1 Half rate constant of methanogens
  
%Global Options
   options = odeset('MaxStep',0.1,'NonNegative',(1:4));
%Calling the Solver
   S = ode45(@groupproject,(0:0.1:tf),c,options);
%add end results
   Rate = [Rate S.y(:,4)];
end
%0.5 times k
for i = 1:12
   K(i) = 0.5;
   %Kinetic Parameters
   Ksa = 20; %mg-S L-1 Half rate constant of anodophilics
   Ksm = 80; %mg-S L-1 Half rate constant of methanogens
%Global Options
   options = odeset('MaxStep',0.1,'NonNegative',(1:4));
%Calling the Solver
   S = ode45(@groupproject,(0:0.1:tf),c,options);
%add end results
   Rate = [Rate S.y(:,4)];
end
%vector subraction for getting percent change
for i = 2:24
   change =[change (((Rate(4*i-3:4*i))-Rate(1:4))./Rate(1:4))];
      
end
change*100;% percent change for all deltas
CH4cath = [2.3e-4, 4.93e-4, 8.5e-4, 1.43e-3, 1.55e-3, 1.54e-3, 0, 7.03e-4, 9.41e-4, 1.38e-03, 0, 2.68e-4, 5.30e-4, 7.37e-4, 2.22e-4, 7.65e-4, 1.25e-3, 1.59e-3, 1.76e-3, 1.95e-3];
ch4t = [1.84, 2.88, 3.98, 5.72, 6.89, 7.91, 1, 2.78, 4.01, 6.91, 1, 2.01, 2.95, 3.74, 1.84, 2.88, 3.98, 5.72, 6.89, 7.91];
t = (1:7);
figure(1)
subplot(131)
xlabel('Time (Days)');
ylabel('Cathode CH4 (mol)');
hold on
plot(t,S.y(3,t))
scatter(ch4t, CH4cath,'or','LineWidth', 1);
legend('Model','Data');
title('Cathode CH4')
hold on
subplot(132)
xlabel('Time (Days)');
ylabel('Anode Acetate (mol)');
hold on
plot(t,S.y(1,t))
subplot(133)
xlabel('Time (Days)');
ylabel('Cathode CO2 (mol)');
hold on
plot(t,S.y(7,t))
figure(2)
bar(change*100)

Answer 1

[Star Strider]

First, do not use global variables! Pass extra parameters as additional parameters. See Passing Extra Parameters for details.

Second, in these assignemtns —

%Anode soluble components
dcdt(1)=(-qa*xa)-(qm*.c(3))+(D*(So-.c(1)));  %Substrate
dcdt(2) = (mua * .c(2)) -1*(Kda*.c(2)) -1*(alpha_a*D*.c(2)); %anodophillic
dcdt(3) = (mum * .c(3)) -1*(Kdm*.c(3)) -1*(alpha_m*D*.c(3)); %methanogenic
%Cathode (Intracellular) Soluble components
dcdt(4) = (-1 * Y * qa) + ((V * (IMFC/(m * F))) * (1/(Vol* .c(2))));

use .* (not *. or -.) because they make no sense. Addition and subtraction are element-wise by definition, so adding dot-operators to those operations is inappropriate.

Third, ‘xa’ and ‘So’ are never defined in ‘rungroupproject’, nor are they included in the global declarations in it (that should not be there anyway), so they are empty and that stops the integration because it is not possible to assign an empty element to a subscripted variable such as ‘dcdt(1)’ (unless it is a cell array element, however that does not apply here).

There could very well be other problems, however I stopped looking after discovering these.