Error using odearguments (line 21) When the first argument to ode23s is a function handle, the tspan argument must have at least two elements.

Hi everyone,

I am trying to solve a set of odes using ode23s, yet I get the error message that states:

Error using odearguments (line 21)
When the first argument to ode23s is a function handle, the tspan argument must have at least two elements.

Error in ode23s (line 121)
    = odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in funcKinetics1 (line 3)
[~, y] = ode23s(@(t,y) funcODE(t, y, p), t, [cs1c0,cs2c0,cs3c0,cs4c0,cs5c0,cs6c0,ce1c0,ce2c0,ce3c0,ce4c0]);

Error in RunFile1 (line 76)
[cs1,cs2,cs3,cs4,cs5,cs6,ce1,ce2,ce3,ce4,pH] = funcKinetics1(t, cs1c0,cs2c0,cs3c0,cs4c0,cs5c0,cs6c0,ce1c0,ce2c0,ce3c0,ce4c0, p);
 
>> 

I am not sure which two elements are reffered to in this error message.

my code is:

function [cs1,cs2,cs3,cs4,cs5,cs6,ce1,ce2,ce3,ce4,pH] = funcKinetics1(t, cs1c0,cs2c0,cs3c0,cs4c0,cs5c0,cs6c0,ce1c0,ce2c0,ce3c0,ce4c0, p)

[~, y] = ode23s(@(t,y) funcODE(t, y, p), t, [cs1c0,cs2c0,cs3c0,cs4c0,cs5c0,cs6c0,ce1c0,ce2c0,ce3c0,ce4c0]);

cs1 = y(:,1);       % High Reactive Lignin
cs2 = y(:,2);       % Low Reactive Lignin
cs3 = y(:,3);       % Cellulose
cs4 = y(:,4);       % Xylan
cs5 = y(:,5);       % Glucomannan
cs6 = y(:,6);       % Acetyls
ce1 = y(:,7);       % Total dissolved lignin
ce2 = y(:,8);       % Total dissolved carbohydrates
ce3 = y(:,9);       % Total acetic acid
ce4 = y(:,10);      % Total sulphite

pH = zeros(length(t), 1);
for i = 1 : length(t)
    [~, pH(i)] = funcOtherVariables(cs1(i),cs2(i),cs3(i),cs4(i),cs5(i),cs6(i),ce1(i),ce2(i),ce3(i),ce4(i), p);
end

end

function dy_dt = funcODE(t, y, p)

cs1 = y(1);       % High Reactive Lignin
cs2 = y(2);       % Low Reactive Lignin
cs3 = y(3);       % Cellulose
cs4 = y(4);       % Xylan
cs5 = y(5);       % Glucomannan
cs6 = y(6);       % Acetyls
ce1 = y(7);       % Total dissolved lignin
ce2 = y(8);       % Total dissolved carbohydrates
ce3 = y(9);       % Total acetic acid
ce4 = y(10);      % Total sulphite

% Calculate the reaction rate.
r = funcOtherVariables(cs1,cs2,cs3,cs4,cs5,cs6,ce1,ce2,ce3,ce4, p);

% Mass balance (derivative terms)
dcs1_dt = r(1,1);
dcs2_dt = r(2,1);
dcs3_dt = r(3,1);
dcs4_dt = r(4,1);
dcs5_dt = r(5,1);
dcs6_dt = r(6,1);
dce1_dt = - p.ve1R1*r(1,1) - p.ve1R2*r(1,1) - p.ve1R3*r(2,1) - p.ve1R4*r(2,1);
dce2_dt = - p.ve2R5*r(3,1) - p.ve2R6*r(4,1) - p.ve2R7*r(5,1);
dce3_dt = - p.ve3R8*r(6,1);
dce4_dt = p.ve4R1*r(1,1) + p.ve4R2*r(1,1) + p.ve4R3*r(2,1) + p.ve4R4*r(2,1)...  
          + p.ve4R5*r(3,1) + p.ve4R6*r(4,1) + p.ve4R7*r(5,1);

% Converting to vector form:
dy_dt = [dcs1_dt;dcs2_dt;dcs3_dt;dcs4_dt;dcs5_dt;dcs6_dt;dce1_dt;dce2_dt;dce3_dt;dce4_dt];

end

function [r, pH] = funcOtherVariables(cs1,cs2,cs3,cs4,cs5,cs6,ce1,ce2,ce3,ce4, p)


pH = fzero(@(pH) funcCharge(pH, cs1,cs2,cs3,cs4,cs5,cs6,ce1,ce2,ce3,ce4, p), 7);

cH = 10^-pH;

% Rate Equations
r1 = -p.k1*cs1^p.vs1R1 * (ce4*((ce4*p.KH2SO3*p.KHSO3)/(cH^2 + cH*p.KH2SO3 + p.KH2SO3*p.KHSO3)))^p.ve4R1 ... 
     -p.k2*cs1^p.vs1R2 * (ce4*((ce4*cH*p.KHSO3)/(cH^2 + cH*p.KH2SO3 + p.KH2SO3*p.KHSO3)))^p.ve4R2;
                       
r2 = -p.k3*cs2^p.vs2R3 * (ce4*((ce4*p.KH2SO3*p.KHSO3)/(cH^2 + cH*p.KH2SO3 + p.KH2SO3*p.KHSO3)))^p.ve4R3 ... 
     -p.k4*cs2^p.vs2R4 * (ce4*((ce4*cH*p.KHSO3)/(cH^2 + cH*p.KH2SO3 + p.KH2SO3*p.KHSO3)))^p.ve4R4;
                       
r3 = -p.k5*cs3^p.vs3R5 * (ce4*((ce4*p.KH2SO3*p.KHSO3)/(cH^2 + cH*p.KH2SO3 + p.KH2SO3*p.KHSO3)))^p.ve4R5; 

r4 = -p.k6*cs4^p.vs4R6 * (ce4*((ce4*cH*p.KHSO3)/(cH^2 + cH*p.KH2SO3 + p.KH2SO3*p.KHSO3)))^p.ve4R6;
                       
r5 = -p.k7*cs5^p.vs5R7 * (ce4*((ce4*p.KH2SO3*p.KHSO3)/(cH^2 + cH*p.KH2SO3 + p.KH2SO3*p.KHSO3)))^p.ve4R7; 

r6 = -p.k8*cs6^p.vs6R8 * cH^p.vHR8 - p.k9*cs6^p.vs6R9 * (p.Kw/cH)^p.vOHR9;

r = [r1;r2;r3;r4;r5;r6];

end

function z = funcCharge(pH, p)

cH = 10^-pH;

% Left-hand side of the charge balance
LHS = p.ce6 + cH;

% Right-hand side of the charge balance
RHS = p.Kw/cH +  2*(s.ce4*p.KH2SO3*p.KHSO3/(cH^2 + cH*p.KH2SO3 + p.KH2SO3*P.KHSO3)) +...
     (s.ce4*cH*p.KH2SO3/(cH^2 + cH*p.KH2SO3 + p.KH2SO3*p.KHSO3))+ ...
     2*(p.ce6*p.KH2CO3*p.KHCO3/(cH^2 + cH*p.KH2CO3 + p.KH2CO3*P.KHCO3)) +...
     (s.ce4*cH*p.KH2CO3/(cH^2 + cH*p.KH2CO3 + p.KH2CO3*p.KHCO3))+ ...
     (s.ce3*p.KAH)/(p.KAH + cH);
 
z = LHS - RHS;  % Should be equal to zero

end

And I am running it using:

%% Parameters

% Rate constants (various units)
p.k1      = 2.95e-1;
p.k2      = 2.0e-2;
p.k3      = 1.0e-1;
p.k4      = 3.0e-1;
p.k5      = 3.0e-3;
p.k6      = 2.2;
p.k7      = 1.0e-2;
p.k8      = 5.0e-2;
p.k9      = 3.0e3;

% Equilibrium constants (various units)
p.KH2CO3    = 1e-3;
p.KHCO3     = 1e-7;
p.KH2SO3    = 1e-3;
p.KHSO3     = 1e-7;
p.KCH2OOH   = 1e-4;
p.KCH3COOH  = 1e-3;
p.Kw        = 1e-8;
p.KAH       = 1e-3;

% Constant concentration species (mol/L)

p.ce5 = 1;       % Total carbonates
p.ce6 = 1;       % Total sodium ions

% Stoichiometric coefficients

p.vs1R1 = 1;
p.vs1R2 = 1;
p.vs2R3 = 1;
p.vs2R4 = 1;
p.vs3R5 = 1;
p.vs4R6 = 1;
p.vs5R7 = 1;
p.vs6R8 = 1;
p.vs6R9 = 1;
p.ve1R1 = 1;
p.ve1R2 = 1;
p.ve1R3 = 1;
p.ve1R4 = 1;
p.ve2R5 = 1;
p.ve2R6 = 1;
p.ve2R7 = 1;
p.ve3R8 = 1;
p.ve4R1 = 1;
p.ve4R2 = 1;
p.ve4R3 = 1;
p.ve4R4 = 1;  
p.ve4R5 = 1;
p.ve4R6 = 1;
p.ve4R7 = 1;
p.vHR8  = 1;
p.vOHR9 = 1;


%--------------------------------------------------------------------------

%% Initial conditions
cs1c0 = 1;       % High Reactive Lignin
cs2c0 = 1;       % Low Reactive Lignin
cs3c0 = 1;       % Cellulose
cs4c0 = 1;       % Xylan
cs5c0 = 1;       % Glucomannan
cs6c0 = 1;       % Acetyls
ce1c0 = 1;       % Total dissolved lignin
ce2c0 = 1;       % Total dissolved carbohydrates
ce3c0 = 1;       % Total acetic acid
ce4c0 = 1;       % Total sulphite
%------------------------------------------------------

%% Numerical integration
t = linspace(0, 100, 1);
[cs1,cs2,cs3,cs4,cs5,cs6,ce1,ce2,ce3,ce4,pH] = funcKinetics1(t, cs1c0,cs2c0,cs3c0,cs4c0,cs5c0,cs6c0,ce1c0,ce2c0,ce3c0,ce4c0, p);

I thank you.

Answer 1

[Bjorn Gustavsson]

Your t variable will only have one component as you define it. My guess is that you want 100 elements between 0 and 1. If that's right you simply swapped the arguments to linspace:

t = linspace(0, 1, 100);

HTH

[Dursman Mchabe]

Oops! Typing error.

[Bjorn Gustavsson]

Yeah, happens to us all - well to me all the time...

...as my father says: first things that go blind on you are the eyes.