I am trying to run this code for parameter estimation, but unable to fix the error. please could you help me to fix the error?

Tem_VL

y0 = 1×13

8381 2000 200 814 800 8029 700 100 14291400 500 150 5000 950

Index exceeds the number of array elements. Index must not exceed 13.

Error in solution>Tem_VL/model_1 (line 108)
+y(3)+y(4)+y(5)+y(6)+y(7)+y(8)+y(9)+y(10))+(c4*b3*y(14))/(y(13)+y(14)))*y(11);

Error in solution>@(t,y)(model_1(t,y,k)) (line 131)
[T Y] = ode23s(@(t,y)(model_1(t,y,k)),tforward,y0);

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

Error in ode23s (line 122)
= odearguments(odeIsFuncHandle, odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error in solution>Tem_VL (line 131)
[T Y] = ode23s(@(t,y)(model_1(t,y,k)),tforward,y0);

function Tem_VL
% This function fits the first set of BSfludata to and SIR model

clear all
close all
clc

% loading VL data 
load VL.txt 

% specifying higher precision
format short

tdata = VL(:,1); % define array with t-coordinates of the data
qdata = VL(:,2); % define array with y-coordinates of the data

tforward = 3:0.01:14; % t mesh for the solution of the differential equation

tmeasure = [1:100:1101]; % selects the points in the solution
% corresponding to the t values of tdata
y0=[8381 2000 200 814 800 8029 700 100 14291400 500 150 5000 950]
% initial values of parameters to be fitted
b1 = 0.56;
b2 = 0.75;
b3 = 0.70;
b4 = 0.56;
c1 = 4;
c2 = 8;
c3 = 0.142;
c4 = 0.333;
d1 = 0.018;
d2 = 0.018;
e = 0.33;
f = 0.33;
h1 = 0.58;
h2 = 0.62;
n1 = 6;
n2 = 6;
n3 = 7.5;
n4 = 8.33;
t1 = 0.5294;
t2 = 0.5294;
u = 0.0000367;
v = 0.00714;
m = 0.16;
w1 = 0.104;
w2 = 0.104;
x = 0.000228;
r = 0.16;
z = 0.8;
p=0.95;
j=0.65;

    function dy = model_1(t,y,k) % DE
        % Assignes the parameters in the DE the current value of the parameters
        b1 = k(1);
        b2 = k(2);
        b3 = k(3);
        b4 = k(4);
        c1 = k(5);
        c2 = k(6);
        c3 = k(7);
        c4 = k(8);
        d1 = k(9);
        d2 = k(10);
        e = k(11);
        f = k(12);
        h1 = k(13);
        h2 = k(14);
        n1 = k(15);
        n2 = k(16);
        n3 = k(17);
        n4 = k(18);
        t1 = k(19);
        t2 = k(20);
        u = k(21);
        v = k(22);
        m = k(23);
        w1 = k(24);
        w2 = k(25);
        x = k(26);
        r = k(27);
        z = k(28);
        p = k(29);
        j = k(30);
        
        
        % assigns zeros to dy 
        dy = zeros(14,1); 
        %NH=y(1)+y(2)+y(3)+y(4)+y(5)+y(6)+y(7)+y(8)+y(9)+y(10)
        dy(1)=n1-(u+(b1*y(12))/(y(1)+y(2)+y(3)+y(4)+y(5)+y(6)+y(7)+y(8)...
            +y(9)+y(10)))*y(1);% RHS of first equation
        dy(2)=((c1*b1*y(12))/(y(1)+y(2)+y(3)+y(4)+y(5)+y(6)+y(7)+y(8)+y(9)...
            +y(10)))*y(1)-(u+r)*y(2);% RHS of second equation
        dy(3)=p*r*y(2)-(e+w1+u)*y(3);
        dy(4)=r*(1-p)*y(2)+e*y(3)-(t1+d1+u)*y(4);
        dy(5)=w1*y(3)+t1*y(4)-u*y(5);
        dy(6)=n2-(u+(c1*b1*y(12))/(y(1)+y(2)+y(3)+y(4)+y(5)+y(6)+y(7)+y(8)...
            +y(9)+y(10)))*y(6);
        dy(7)=((c1*b1*y(12))/(y(1)+y(2)+y(3)+y(4)+y(5)+y(6)+y(7)+y(8)+y(9)...
            +y(10)))*y(6)-(u+m)*y(7);
        dy(8)=m*j*y(7)-(f+w2+u)*y(8);
        dy(9)=m*(1-j)*y(7)+f*y(8)-(t2+d2+u)*y(9);
        dy(10)=w2*y(8)+t2*y(9)-u*y(10);
        dy(11)=n3-(v+c3*b2*(y(3)+h1*y(8)+z*(y(4)+h2*y(9)))/(y(1)+y(2)...
            +y(3)+y(4)+y(5)+y(6)+y(7)+y(8)+y(9)+y(10))+(c4*b3*y(14))/(y(13)+y(14)))*y(11);
        dy(12)=(c3*b2*(y(3)+h1*y(8)+z*(y(4)+h2*y(9)))/(y(1)+y(2)+y(3)+y(4)...
            +y(5)+y(6)+y(7)+y(8)+y(9)+y(10))+(c4*b3*y(14))/(y(13)...
            +y(14)))*y(11)-v*y(12);
        dy(13)=n4-(x+(c2*b4*y(12))/(y(13)+y(14)))*y(13);
        dy(14)=((c2*b4*y(12))/(y(13)+y(14)))*y(13)-x*y(14);
    end

%computing the error in the data
    function error_in_data = moder(k) 
        [T Y] = ode23s(@(t,y)(model_1(t,y,k)),tforward,y0);
        % solves the DE; output is written in T and Y
        % assignts the y-coordinates of the solution at the t-coordinates of tdata
        q = Y(tmeasure(:),2); 
        
        %computes SSE
        error_in_data = sum((q - qdata).^2) 
        
    end
% main routine; assigns initial values of parameters
k = [b1 b2 b3 b4 c1 c2 c3 c4 d1 d2 e f h1 h2 n1 n2 n3 n4 ...
    t1 t2 u v m w1 w2 x r z p j];

[T Y] = ode23s(@(t,y)(model_1(t,y,k)),tforward,y0);

% solves the DE with the initial values of the parameters

yint = Y(tmeasure(:),2);

% assigns the y-coordinates of the solution at tdata to yint

figure(1)
subplot(1,2,1);
plot(tdata,qdata,'r.');
hold on
plot(tdata,yint,'g-'); % plotting of solution and data with initial
% guesses for the parameters
xlabel('time in days');
ylabel('Number of cases');
axis([3 14 0 350]);
[k,fval] = fminsearch(@moder,k); % minimization routine; assigns the new
% values of parameters to k and the SSE
% to fval

disp(k);
[T Y] = ode23s(@(t,y)(model_1(t,y,k)),tforward,y0);
% solving the DE with the final values of the
% parameters

yint = Y(tmeasure(:),2); % computing the y-coordinates corresponding to the
% tdata

subplot(1,2,2)
plot(tdata,qdata,'r.');
hold on
plot(tdata,yint,'c-');
xlabel('time in days'); % plotting final fit
ylabel('Number of cases');
axis([3 20 0 350]);
end

[Dyuman Joshi]

Please format your code properly, and state the full error message.

[Temesgen]

@Dyuman Joshi the following are the error i have found

Attempted to access y(14); index out of bounds because numel(y)=13.

Error in Tem_VL/model_1 (line 103)

dy(11)=n3-(v+c3*b2*(y(3)+h1*y(8)+z*(y(4)+h2*y(9)))/(y(1)+y(2)...

Error in Tem_VL>@(t,y)(model_1(t,y,k)) (line 127)

[T Y] = ode23s(@(t,y)(model_1(t,y,k)),tforward,y0);

Error in odearguments (line 87)

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

Error in ode23s (line 119)

[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in Tem_VL (line 127)

[T Y] = ode23s(@(t,y)(model_1(t,y,k)),tforward,y0);

[Walter Roberson]

We do not have your data to test with.

[Temesgen]

Dear Walter thank you very much for your correction. Here is my data. but I found this error again. I hope you are fixing my problem. I am waiting you!!

Error using plot

Vectors must be the same length.

Error in Tem_VL (line 139)

plot(tdata,yint,'g-'); % plotting of solution and data with initial

Answer 1

[Walter Roberson]

y0=[8381 2000 200 814 800 8029 700 100 14291400 500 150 5000 950]

That is 13 entries, one of which is much larger value than the rest. We might wonder, for example, whether the data is intended to be

y0=[8381 2000 200 814 800 8029 700 100 1429 1400 500 150 5000 950]

??

        dy = zeros(14,1); 
 

Your code is obviously expecting to work with 14 state variables, but you are only passing in 13 state values.

[Walter Roberson]

Tem_VL

y0 = 1×14

8381 2000 200 814 800 8029 700 100 1429 1400 500 150 5000 950

Error using plot
Vectors must be the same length.

Error in solution>Tem_VL (line 143)
plot(tdata,yint,'g-'); % plotting of solution and data with initial

function Tem_VL
% This function fits the first set of BSfludata to and SIR model

clear all
close all
clc

% loading VL data 
load VL.txt 

% specifying higher precision
format short

tdata = VL(:,1); % define array with t-coordinates of the data
qdata = VL(:,2); % define array with y-coordinates of the data

tforward = 3:0.01:14; % t mesh for the solution of the differential equation

tmeasure = [1:100:1101]; % selects the points in the solution
% corresponding to the t values of tdata
y0=[8381 2000 200 814 800 8029 700 100 1429 1400 500 150 5000 950]
% initial values of parameters to be fitted
b1 = 0.56;
b2 = 0.75;
b3 = 0.70;
b4 = 0.56;
c1 = 4;
c2 = 8;
c3 = 0.142;
c4 = 0.333;
d1 = 0.018;
d2 = 0.018;
e = 0.33;
f = 0.33;
h1 = 0.58;
h2 = 0.62;
n1 = 6;
n2 = 6;
n3 = 7.5;
n4 = 8.33;
t1 = 0.5294;
t2 = 0.5294;
u = 0.0000367;
v = 0.00714;
m = 0.16;
w1 = 0.104;
w2 = 0.104;
x = 0.000228;
r = 0.16;
z = 0.8;
p=0.95;
j=0.65;

    function dy = model_1(t,y,k) % DE
        % Assignes the parameters in the DE the current value of the parameters
        b1 = k(1);
        b2 = k(2);
        b3 = k(3);
        b4 = k(4);
        c1 = k(5);
        c2 = k(6);
        c3 = k(7);
        c4 = k(8);
        d1 = k(9);
        d2 = k(10);
        e = k(11);
        f = k(12);
        h1 = k(13);
        h2 = k(14);
        n1 = k(15);
        n2 = k(16);
        n3 = k(17);
        n4 = k(18);
        t1 = k(19);
        t2 = k(20);
        u = k(21);
        v = k(22);
        m = k(23);
        w1 = k(24);
        w2 = k(25);
        x = k(26);
        r = k(27);
        z = k(28);
        p = k(29);
        j = k(30);
        
        
        % assigns zeros to dy 
        dy = zeros(14,1); 
        %NH=y(1)+y(2)+y(3)+y(4)+y(5)+y(6)+y(7)+y(8)+y(9)+y(10)
        dy(1)=n1-(u+(b1*y(12))/(y(1)+y(2)+y(3)+y(4)+y(5)+y(6)+y(7)+y(8)...
            +y(9)+y(10)))*y(1);% RHS of first equation
        dy(2)=((c1*b1*y(12))/(y(1)+y(2)+y(3)+y(4)+y(5)+y(6)+y(7)+y(8)+y(9)...
            +y(10)))*y(1)-(u+r)*y(2);% RHS of second equation
        dy(3)=p*r*y(2)-(e+w1+u)*y(3);
        dy(4)=r*(1-p)*y(2)+e*y(3)-(t1+d1+u)*y(4);
        dy(5)=w1*y(3)+t1*y(4)-u*y(5);
        dy(6)=n2-(u+(c1*b1*y(12))/(y(1)+y(2)+y(3)+y(4)+y(5)+y(6)+y(7)+y(8)...
            +y(9)+y(10)))*y(6);
        dy(7)=((c1*b1*y(12))/(y(1)+y(2)+y(3)+y(4)+y(5)+y(6)+y(7)+y(8)+y(9)...
            +y(10)))*y(6)-(u+m)*y(7);
        dy(8)=m*j*y(7)-(f+w2+u)*y(8);
        dy(9)=m*(1-j)*y(7)+f*y(8)-(t2+d2+u)*y(9);
        dy(10)=w2*y(8)+t2*y(9)-u*y(10);
        dy(11)=n3-(v+c3*b2*(y(3)+h1*y(8)+z*(y(4)+h2*y(9)))/(y(1)+y(2)...
            +y(3)+y(4)+y(5)+y(6)+y(7)+y(8)+y(9)+y(10))+(c4*b3*y(14))/(y(13)+y(14)))*y(11);
        dy(12)=(c3*b2*(y(3)+h1*y(8)+z*(y(4)+h2*y(9)))/(y(1)+y(2)+y(3)+y(4)...
            +y(5)+y(6)+y(7)+y(8)+y(9)+y(10))+(c4*b3*y(14))/(y(13)...
            +y(14)))*y(11)-v*y(12);
        dy(13)=n4-(x+(c2*b4*y(12))/(y(13)+y(14)))*y(13);
        dy(14)=((c2*b4*y(12))/(y(13)+y(14)))*y(13)-x*y(14);
    end

%computing the error in the data
    function error_in_data = moder(k) 
        [T Y] = ode23s(@(t,y)(model_1(t,y,k)),tforward,y0);
        % solves the DE; output is written in T and Y
        % assignts the y-coordinates of the solution at the t-coordinates of tdata
        q = Y(tmeasure(:),2); 
        
        %computes SSE
        error_in_data = sum((q - qdata).^2) 
        
    end
% main routine; assigns initial values of parameters
k = [b1 b2 b3 b4 c1 c2 c3 c4 d1 d2 e f h1 h2 n1 n2 n3 n4 ...
    t1 t2 u v m w1 w2 x r z p j];

[T Y] = ode23s(@(t,y)(model_1(t,y,k)),tforward,y0);

% solves the DE with the initial values of the parameters

yint = Y(tmeasure(:),2);

% assigns the y-coordinates of the solution at tdata to yint

figure(1)
subplot(1,2,1);
plot(tdata,qdata,'r.');
hold on
plot(tdata,yint,'g-'); % plotting of solution and data with initial
% guesses for the parameters
xlabel('time in days');
ylabel('Number of cases');
axis([3 14 0 350]);
[k,fval] = fminsearch(@moder,k); % minimization routine; assigns the new
% values of parameters to k and the SSE
% to fval

disp(k);
[T Y] = ode23s(@(t,y)(model_1(t,y,k)),tforward,y0);
% solving the DE with the final values of the
% parameters

yint = Y(tmeasure(:),2); % computing the y-coordinates corresponding to the
% tdata

subplot(1,2,2)
plot(tdata,qdata,'r.');
hold on
plot(tdata,yint,'c-');
xlabel('time in days'); % plotting final fit
ylabel('Number of cases');
axis([3 20 0 350]);
end

[Walter Roberson]

tdata = VL(:,1); % define array with t-coordinates of the data

That is whatever size is read from the file.

tforward = 3:0.01:14

That is probably 1101 elements.

[T Y] = ode23s(@(t,y)(model_1(t,y,k)),tforward,y0);

You pass those 1101 elements as the times to return for ode23s. Those will be returned into T (unless the ode23s has to end early)

tmeasure = [1:100:1101]; % selects the points in the solution

That is 11 points.

yint = Y(tmeasure(:),2); % computing the y-coordinates corresponding to the

That is going to be 11 points.

plot(tdata,yint,'g-'); % plotting of solution and data with initial

tdata is whatever size was read from the file. yint is 11 points selected out of the 1101 generated for the ode. Those are probably not going to be the same size.

The values in tdata could potentially not intersect 3:0.01:14 at all, or might overlap, or might be fully contained inside (as outside observers we do not know.) You could consider interpolating the values in tdata against the values in tforward to find the closing matching indices and extract the corresponding Y(:,2) values, and plotting tdata against those values.

Or in the case where you are certain that tdata (from the file) is entirely contained within 3:0.01:14 you could create tforward as being [3 tdata] and then throwing away the first row of result (corresponding to the time 3), and the remanining rows would correspond to data. Unless, that is, tdata happens to start with 3 itself. Or unless it happens to start before 3...

[Temesgen]

Dear Walter, thank you very much again for your swift reply with your valuable direction. Oucourse, this program is working for simple sample model and by considering this, I modified that to this complex model. I am a beginer to coding in Matlab and my first time to post on hear.

I would appreciate you if you could observe the sample code attached herewith and tell me the suggested source of error during modification??