WATERMODEL returns a vector of length 825, but the length of initial conditions vector is 3. The vector returned by WATERMODEL and the initial conditions vector must have the

I am trying to solve the ODEs defined by equations:

dL1_dt=-fTr1-fDr1;

dL2_dt=fDr1-fTr2-fDr2;

dL3_dt=fDr2-fTr3-fDr3;

But I am getting an error, I cant figure out what is wrong. Parameters used in calculating the right hand sides of the equations change after every time step, depending on the data contained in "dummyData.csv". I attach the file

tspan=1:275;  % [d] 275 is length of data in data table
%Initial conditions    
L0=[54;80;144];  %[mm], 
[t,L]=ode45(@waterModel,tspan,L0);

Error using odearguments
WATERMODEL returns a vector of length 825, but the length of initial conditions vector is 3. The vector returned by WATERMODEL and the initial conditions vector must have the same number of
elements.

Error in ode45 (line 107)
odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

%% ODE FUNCTION
function  dLdt=waterModel(t,Y)
L1=Y(1);
L2=Y(2);
L3=Y(3);
% Load file with variables temperature, ET etc,
%%
data=readtable("dummyData.csv");

% Use the data to calculate waterparamters 

p=waterParameters(data.Temperature); % Call to the nested function for calculating water parameters
% p is a structure of parameters
%% Plant transpiration
   WAI1=(L1-p.pwp(1))./(p.fc(1)-p.pwp(1));
   WAI2=(L2-p.pwp(2))./(p.fc(2)-p.pwp(2));
   WAI3=(L3-p.pwp(3))./(p.fc(3)-p.pwp(3));
    if WAI1<p.WAIc
         krt1=WAI1./p.WAIc;
    else
         krt1=1;
    end
   
    if WAI2<p.WAIc
            krt2=WAI2./p.WAIc;
    else
            krt2=1;
    end     
    if WAI3<p.WAIc
            krt3=WAI3./p.WAIc;
    else
            krt3=1;
    end
%% Transpiration from each soil layer 1-3[mm]
        fTr1=krt1.*p.kra.*p.krf1.*data.ET;
        fTr2=krt2.*p.kra.*p.krf2.*data.ET;
        fTr3=krt3.*p.kra.*p.krf3.*data.ET;
        
%% Soil water Drainage from both layers
     if L1>p.fc(1)
           fDr1=L1-p.fc(1);
     else
            fDr1=0;
     end    
     if L2>p.fc(2)
            fDr2=L2-p.fc(2);
     else
            fDr2=0;
     end    
     if L3>p.fc(3)
            fDr3=L3-p.fc(3);
     else
            fDr3=0;
     end
%% Total Water availability index
%WAI=(L1-p.pwp(1)+L2-p.pwp(2)+L3-p.pwp(3))./(p.fc(1)-p.pwp(1)+p.fc(2)-p.pwp(2)+p.fc(3)-p.pwp(3));
%WAI = max(WAI,0);  
%%    Controlled inputs
%firr = 0;         % [mm d-1] Irrigation

%% Differential equations    
dL1_dt=-fTr1-fDr1; % dL1_dt=fpe-fTr1-fEv-fDr1
dL2_dt=fDr1-fTr2-fDr2;
dL3_dt=fDr2-fTr3-fDr3;
dLdt=[dL1_dt;dL2_dt;dL3_dt];

function p=waterParameters(Tgrh)
% Function for calculating water parameters, takes temperature of
% greenhouse [Tgrh] and returns a structure of parametrs
p.WAIc=0.75;              %    [-] WAI critical, range (0.5-0.8)
p.theta_fc1=0.36;         %    [-] Field capacity of soil layer 1
p.theta_fc2=0.32;         %    [-] Field capacity of soil layer 2
p.theta_fc3=0.24;         %    [-] Field capacity of soil layer 3
p.theta_pwp1= 0.21;       %    [-] Permanent wilting point of soil layer 1   
p.theta_pwp2=0.17;        %    [-] Permanent wilting point of soil layer 2   
p.theta_pwp3=0.1;         %    [-] Permanent wilting point of soil layer 3
p.D1=150;                 %    [mm] Depth of Soil layer 1
p.D2=250;                 %    [mm] Depth of Soil layer 2
p.D3=600;                 %    [mm] Depth of Soil layer 3
p.krf1=0.25;              %    [-] Rootfraction layer 1 (guess)
p.krf2=0.50;              %    [-] Rootfraction layer 2 (guess)
p.krf3=0.25;              %    [-] Rootfraction layer 3 (guess)
p.kra= 0.0408.* exp(0.19.*Tgrh); % [-] Root activity, varies with Tgrh (so its a vector of kra per time step)
% [mm] Field capacities of both soil layers
p.fc=[p.theta_fc1*p.D1, p.theta_fc2*p.D2, p.theta_fc3*p.D3];

% [mm] Permanent wilting points
p.pwp = [p.theta_pwp1*p.D1, p.theta_pwp2*p.D2, p.theta_pwp3*p.D3]; 
 
end
end

Answer 1

[Torsten]

% Load file with variables temperature, ET etc,
%%
data=readtable("https://de.mathworks.com/matlabcentral/answers/uploaded_files/1150115/dummyData.csv");
tspan=data.TIME(1):data.TIME(end);  % [d] 275 is length of data in data table
%Initial conditions    
L0=[54;80;144];  %[mm], 
[t,L]=ode45(@(t,y)waterModel(t,y,data),tspan,L0);
plot(t,L)

%% ODE FUNCTION
function  dLdt=waterModel(t,Y,data)
L1=Y(1);
L2=Y(2);
L3=Y(3);
% Use the data to calculate waterparamters 

p=waterParameters(data.TIME,data.Temperature,t); % Call to the nested function for calculating water parameters
% p is a structure of parameters
%% Plant transpiration
   WAI1=(L1-p.pwp(1))./(p.fc(1)-p.pwp(1));
   WAI2=(L2-p.pwp(2))./(p.fc(2)-p.pwp(2));
   WAI3=(L3-p.pwp(3))./(p.fc(3)-p.pwp(3));
    if WAI1<p.WAIc
         krt1=WAI1./p.WAIc;
    else
         krt1=1;
    end
   
    if WAI2<p.WAIc
            krt2=WAI2./p.WAIc;
    else
            krt2=1;
    end     
    if WAI3<p.WAIc
            krt3=WAI3./p.WAIc;
    else
            krt3=1;
    end
%% Transpiration from each soil layer 1-3[mm]
        fTr1=krt1.*p.kra.*p.krf1.*interp1(data.TIME,data.ET,t);
        fTr2=krt2.*p.kra.*p.krf2.*interp1(data.TIME,data.ET,t);
        fTr3=krt3.*p.kra.*p.krf3.*interp1(data.TIME,data.ET,t);
        
%% Soil water Drainage from both layers
     if L1>p.fc(1)
           fDr1=L1-p.fc(1);
     else
            fDr1=0;
     end    
     if L2>p.fc(2)
            fDr2=L2-p.fc(2);
     else
            fDr2=0;
     end    
     if L3>p.fc(3)
            fDr3=L3-p.fc(3);
     else
            fDr3=0;
     end
%% Total Water availability index
%WAI=(L1-p.pwp(1)+L2-p.pwp(2)+L3-p.pwp(3))./(p.fc(1)-p.pwp(1)+p.fc(2)-p.pwp(2)+p.fc(3)-p.pwp(3));
%WAI = max(WAI,0);  
%%    Controlled inputs
%firr = 0;         % [mm d-1] Irrigation

%% Differential equations    
dL1_dt=-fTr1-fDr1; % dL1_dt=fpe-fTr1-fEv-fDr1
dL2_dt=fDr1-fTr2-fDr2;
dL3_dt=fDr2-fTr3-fDr3;
dLdt=[dL1_dt;dL2_dt;dL3_dt];

function p=waterParameters(Time,Tgrh,t)
% Function for calculating water parameters, takes temperature of
% greenhouse [Tgrh] and returns a structure of parametrs
p.WAIc=0.75;              %    [-] WAI critical, range (0.5-0.8)
p.theta_fc1=0.36;         %    [-] Field capacity of soil layer 1
p.theta_fc2=0.32;         %    [-] Field capacity of soil layer 2
p.theta_fc3=0.24;         %    [-] Field capacity of soil layer 3
p.theta_pwp1= 0.21;       %    [-] Permanent wilting point of soil layer 1   
p.theta_pwp2=0.17;        %    [-] Permanent wilting point of soil layer 2   
p.theta_pwp3=0.1;         %    [-] Permanent wilting point of soil layer 3
p.D1=150;                 %    [mm] Depth of Soil layer 1
p.D2=250;                 %    [mm] Depth of Soil layer 2
p.D3=600;                 %    [mm] Depth of Soil layer 3
p.krf1=0.25;              %    [-] Rootfraction layer 1 (guess)
p.krf2=0.50;              %    [-] Rootfraction layer 2 (guess)
p.krf3=0.25;              %    [-] Rootfraction layer 3 (guess)
p.kra= 0.0408.* exp(0.19.*interp1(Time,Tgrh,t)); % [-] Root activity, varies with Tgrh (so its a vector of kra per time step)
% [mm] Field capacities of both soil layers
p.fc=[p.theta_fc1*p.D1, p.theta_fc2*p.D2, p.theta_fc3*p.D3];

% [mm] Permanent wilting points
p.pwp = [p.theta_pwp1*p.D1, p.theta_pwp2*p.D2, p.theta_pwp3*p.D3]; 
 
end
end

[Nick Austinos]

Thank you very much @Torsten it works perfectly Now; I did not no this trick................

[Torsten]

For any time t, you must calculate the suitable values from your data according to the time of integration.

You can't supply the complete column of the data vector.

And read the data once at the beginning of your code and pass them to the ODE function as done in the code above. It would be a waste of time to read them again every time the ODE function is called.