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

function run_LE_FO_a(ne,ext_fcn,t_start,h_norm,t_end,x_start,h,q,a_min,a_max,n);
figure();
hold on;
ne=4;
ext_fcn=@LE_RF_a;
t_start=0;h_norm=0.02;t_end=1;
x_start=[0 0 1 1];h=0.2;q=1;
a_max=-11;
a_min=-12.2;
n=800;
p_step=(a_max-a_min)/n
a=a_min;
% for  p1=1.1:0.1:1.3;
while a<=a_max
    lp=FO_Lyapunov_a(ne,ext_fcn,0,0.2,1,[0 0 1 1],0.2,1,a)
    a=a+p_step;
    plot(a,lp,'.')
    drawnow();
    % end
end
function LE=FO_Lyapunov_a(ne,ext_fcn,t_start,h_norm,t_end,x_start,h,q,a);
% Memory allocation
x=zeros(ne*(ne+1),1); 
x0=x;
% y0=x_start
c=zeros(ne,1);
gsc=c; zn=c;
n_it = round((t_end-t_start)/h_norm);
% Initial values
x(1:ne)=x_start;
i=1;
while i<=ne
    x((ne+1)*i)=1.0;
    i=i+1;
end
t=t_start;
% Main loop
it=1;
while it<=n_it
    % Solutuion of extended ODE system of FO using FDE12 routine
    [T,Y] = FDE12(q,ext_fcn,t,t+h_norm,[0 0 1 1],h,a); 
    t=t+h_norm;
    Y=transpose(Y);
    x=Y(size(Y,1),:); %solution at t+h_norm
    i=1;
    while i<=ne
        j=1;
        while j<=ne;
            x0(ne*i+j)=x(ne*j+i);
            j=j+1;
        end;
        i=i+1;
    end;
    %       construct new orthonormal basis by gram-schmidt
    zn(1)=0.0;
    j=1;
    while j<=ne
        zn(1)=zn(1)+x0(ne*j+1)^2; 
        j=j+1;
    end;
    zn(1)=sqrt(zn(1));
    j=1;
    while j<=ne
        x0(ne*j+1)=x0(ne*j+1)/zn(1); 
        j=j+1;
    end
    j=2;
    while j<=ne
        k=1;
        while k<=j-1
            gsc(k)=0.0;
            l=1;
            while l<=ne;
                gsc(k)=gsc(k)+x0(ne*l+j)*x0(ne*l+k);
                l=l+1;
            end
            k=k+1;
        end
        k=1;
        while k<=ne
            l=1;
            while l<=j-1
                x0(ne*k+j)=x0(ne*k+j)-gsc(l)*x0(ne*k+l);
                l=l+1;
            end
            k=k+1;
        end;
        zn(j)=0.0;
        k=1;
        while k<=ne
            zn(j)=zn(j)+x0(ne*k+j)^2;
            k=k+1;
        end
        zn(j)=sqrt(zn(j));
        k=1;
        while k<=ne
            x0(ne*k+j)=x0(ne*k+j)/zn(j);
            k=k+1;
        end
        j=j+1;
    end
    %       update running vector magnitudes
    k=1;
    while k<=ne;
        c(k)=c(k)+log(zn(k));
        k=k+1;
    end;
    %       normalize exponent        
    k=1;
    while k<=ne
        LE(k)=c(k)/(t-t_start);
        k=k+1;
    end
    i=1;
    while i<=ne
           j=1;
            while j<=ne;
                x(ne*j+i)=x0(ne*i+j);
                j=j+1;
            end
            i=i+1;
        end;
        x=transpose(x);
        it=it+1;
end
function f=LE_RF_a(t,x,a)
%p is the parameter
f=zeros(16,1);
b=0.4;c=11;d=6;e=130;f1=10;
x=X(1); y=X(2); z=X(3);w=X(4);
Y= [X(5), X(9), X(13),  X(17);
    X(6), X(10), X(14), X(18);
    X(7), X(11), X(15), X(19);
    X(8), X(12),  X(16), X(20) ];
f(1)=a*y+(0.2+0.2*(abs(w))).*z;
f(2)=b*((w.^2)-13).*z-c*y;
f(3)=-d*x-e*y-f1*z;
f(4)=(z.^2)-w.^2;%Linearized system
 Jac=[0 a 0.2+0.2*abs(w) 0.2*w.*z/abs(w);
     0 -c b*((w.^2)-13) b*((w*2)-13).*z;
     -d -e -f1 0;
     0 0 2*z -2*w];
  
f(5:20)=Jac*Y;

error:>> run_LE_FO_a

p_step =

1.499999999999999e-03

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

Error in LE_RF_p1 (line 4)

X= [x(4) x(7) x(10);

Error in FDE12>f_vectorfield (line 300)

f = feval(Probl.fdefun,t,y,Probl.param) ;

Error in FDE12 (line 114)

f_temp = f_vectorfield(t0,y0(:,1),Probl) ;

Error in FO_Lyapunov_a (line 21)

[T,Y] = FDE12(q,ext_fcn,t,t+h_norm,[0 0 1 1],h,a);

Error in run_LE_FO_a (line 16)

lp=FO_Lyapunov_a(ne,ext_fcn,0,0.2,1,[0 0 1 1],0.2,1,a)

>>

Answer 1

[Cris LaPierre]

See the full error message:

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

Error in LE_RF_p1 (line 4)

X= [x(4) x(7) x(10);

Apparently x is a scalar, not a vector, meaning it only has a single value. There error is because your code is requesting the 4th, 7th and 10th values, which don't exist.

x = 3;
% works
x(1)

ans = 3

% Your error
x(4)

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

[nune pratyusha]

i didn't get it and i am getting error

>> run_LE_FO_a

p_step =

1.499999999999999e-03

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

Error in LE_RF_a (line 5)

x=X(1); y=X(2); z=X(3);w=X(4);

Error in FDE12>f_vectorfield (line 300)

f = feval(Probl.fdefun,t,y,Probl.param) ;

Error in FDE12 (line 114)

f_temp = f_vectorfield(t0,y0(:,1),Probl) ;

Error in FO_Lyapunov_a (line 21)

[T,Y] = FDE12(q,ext_fcn,t,t+h_norm,[0 0 1 1],h,a);

Error in run_LE_FO_a (line 16)

lp=FO_Lyapunov_a(ne,ext_fcn,0,0.2,1,[0 0 1 1],0.2,1,a)

where i have to modify i didn't understand

[Torsten]

There does not exist an array X in the function LE_RF_a.

Thus setting

x=X(1); y=X(2); z=X(3);w=X(4);

gives an error.

[nune pratyusha]

I changed x,y & z as follows after that also showing same error

function f=LE_RF_a(t,x,a)
%p is the parameter
% x=1:20;
f=zeros(16,1);
b=0.4;c=11;d=6;e=130;f1=10;
X= [x(5) x(9) x(13)  x(17);
    x(6) x(10) x(14) x(18);
    x(7) x(11) x(15) x(19);
    x(8) x(12)  x(16) x(20) ];
f(1)=a*x(2)+(0.2+0.2*(abs(x(4)))).*x(3);
f(2)=b*((x(4).^2)-13).*x(3)-c*x(2);
f(3)=-d*x(1)-e*x(2)-f1*x(3);
f(4)=(x(3).^2)-x(4).^2;
Jac=[0 a 0.2+0.2*abs(x(4)) 0.2*x(4).*x(3)/abs(x(3));
    0 -c b*((x(4).^2)-13) b*((x(4)*2)-13).*x(3);
    -d -e -f1 0;
    0 0 2*x(3) -2*x(4)];

f(5:20)=Jac*X;

error:

>> run_LE_FO_a

p_step =

1.499999999999999e-03

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

Error in LE_RF_a (line 6)

X= [x(5) x(9) x(13) x(17);

Error in FDE12>f_vectorfield (line 300)

f = feval(Probl.fdefun,t,y,Probl.param) ;

Error in FDE12 (line 114)

f_temp = f_vectorfield(t0,y0(:,1),Probl) ;

Error in FO_Lyapunov_a (line 23)

[T,Y] = FDE12(q,ext_fcn,t,t+h_norm,[0 0 1 1],h,a);

Error in run_LE_FO_a (line 16)

lp=FO_Lyapunov_a(4,ext_fcn,0,0.2,1,[0 0 1 1],0.2,1,a)

>>

[Cris LaPierre]

You need to share LE_RF_p1 as that is where the error is occuring. It would be helpful is you also shared all dependent code so that we can run it and recreate the error message.

[nune pratyusha]

function f=LE_RF_a(t,x,a)
%p is the parameter
% x=1:20;
f=zeros(16,1);
b=0.4;c=11;d=6;e=130;f1=10;
X= [x(5) x(9) x(13)  x(17);
    x(6) x(10) x(14) x(18);
    x(7) x(11) x(15) x(19);
    x(8) x(12)  x(16) x(20) ];
f(1)=a*x(2)+(0.2+0.2*(abs(x(4)))).*x(3);
f(2)=b*((x(4).^2)-13).*x(3)-c*x(2);
f(3)=-d*x(1)-e*x(2)-f1*x(3);
f(4)=(x(3).^2)-x(4).^2;
Jac=[0 a 0.2+0.2*abs(x(4)) 0.2*x(4).*x(3)/abs(x(3));
    0 -c b*((x(4).^2)-13) b*((x(4)*2)-13).*x(3);
    -d -e -f1 0;
    0 0 2*x(3) -2*x(4)];

f(5:20)=Jac*X;
function LE=FO_Lyapunov_a(ne,ext_fcn,t_start,h_norm,t_end,x_start,h,q,a);
% Memory allocation
ne=4;
x=zeros(ne*(ne+1),1); 
x0=x;
q=1;
% y0=x_start
c=zeros(ne,1);
gsc=c; zn=c;
n_it = round((t_end-t_start)/h_norm);
% Initial values
x(1:ne)=x_start;
i=1;
while i<=ne
    x((ne+1)*i)=1.0;
    i=i+1;
end
t=t_start;
% Main loop
it=1;
while it<=n_it
    % Solutuion of extended ODE system of FO using FDE12 routine
    [T,Y] = FDE12(q,ext_fcn,t,t+h_norm,[0 0 1 1],h,a); 
    t=t+h_norm;
    Y=transpose(Y);
    x=Y(size(Y,1),:); %solution at t+h_norm
    i=1;
    while i<=ne
        j=1;
        while j<=ne;
            x0(ne*i+j)=x(ne*j+i);
            j=j+1;
        end;
        i=i+1;
    end;
    %       construct new orthonormal basis by gram-schmidt
    zn(1)=0.0;
    j=1;
    while j<=ne
        zn(1)=zn(1)+x0(ne*j+1)^2; 
        j=j+1;
    end;
    zn(1)=sqrt(zn(1));
    j=1;
    while j<=ne
        x0(ne*j+1)=x0(ne*j+1)/zn(1); 
        j=j+1;
    end
    j=2;
    while j<=ne
        k=1;
        while k<=j-1
            gsc(k)=0.0;
            l=1;
            while l<=ne;
                gsc(k)=gsc(k)+x0(ne*l+j)*x0(ne*l+k);
                l=l+1;
            end
            k=k+1;
        end
        k=1;
        while k<=ne
            l=1;
            while l<=j-1
                x0(ne*k+j)=x0(ne*k+j)-gsc(l)*x0(ne*k+l);
                l=l+1;
            end
            k=k+1;
        end;
        zn(j)=0.0;
        k=1;
        while k<=ne
            zn(j)=zn(j)+x0(ne*k+j)^2;
            k=k+1;
        end
        zn(j)=sqrt(zn(j));
        k=1;
        while k<=ne
            x0(ne*k+j)=x0(ne*k+j)/zn(j);
            k=k+1;
        end
        j=j+1;
    end
    %       update running vector magnitudes
    k=1;
    while k<=ne;
        c(k)=c(k)+log(zn(k));
        k=k+1;
    end;
    %       normalize exponent        
    k=1;
    while k<=ne
        LE(k)=c(k)/(t-t_start);
        k=k+1;
    end
    i=1;
    while i<=ne
        j=1;
        while j<=ne;
            x(ne*j+i)=x0(ne*i+j);
            j=j+1;
        end
        i=i+1;
    end;
    x=transpose(x);
    it=it+1;
end
function run_LE_FO_a(ne,ext_fcn,t_start,h_norm,t_end,x_start,h,q,a_min,a_max,n);
figure();
hold on;
ne=4;
ext_fcn=@LE_RF_a;
t_start=0;h_norm=0.02;t_end=1;
x_start=[0 0 1 1];h=0.2;q=1;

a_max=-11;
a_min=-12.2;
n=800;
p_step=(a_max-a_min)/n
a=a_min;
% for  p1=1.1:0.1:1.3;
while a<=a_max
    lp=FO_Lyapunov_a(4,ext_fcn,0,0.2,1,[0 0 1 1],0.2,1,a)
    a=a+p_step;
    plot(a,lp,'.')
    drawnow();
    % end
end

[Cris LaPierre]

LE_RF_p1 is still missing.

Also, how do you call FDE12? There are a lot of inputs we do not know.

FDE12(alpha,fdefun,t0,tfinal,y0,h,param,mu,mu_tol)

[nune pratyusha]

[T,Y] = FDE12(q,ext_fcn,t,t+h_norm,[0 0 1 1],h,a);

using above command i call FDE12

and in the place of LE_RF_p1 i changed LE_RF_a

[Cris LaPierre]

The x input to LE_RF_a is a single value, not a vector. The error is happening in the following line of code in FDE12

% Check compatibility size of the problem with size of the vector field

f_temp = f_vectorfield(t0,y0(:,1),Probl) ;

This in turn calls the following function

f = feval(Probl.fdefun,t,y,Probl.param)

Here, y corresponds to the input 'x' in f=LE_RF_a(t,x,a). You are passing in a single value (y0(:,1) where y0 = [0 0 1 1]). Even if you didn't index y0, it only has 4 values. Your code in LE_RF_a is written assuming x has at least 20 elements.

The values of y0 are hard coded into your call to FDE12 inside FO_Lyapunov_a:

% Solutuion of extended ODE system of FO using FDE12 routine

[T,Y] = FDE12(q,ext_fcn,t,t+h_norm,[0 0 1 1],h,a);

Here's an executable version of your code with the error for others who may want to help.

ne=4;
ext_fcn=@LE_RF_a;
t_start=0;h_norm=0.02;t_end=1;
x_start=[0 0 1 1];h=0.2;q=1;
a_max=-11;
a_min=-12.2;
n=800;
p_step=(a_max-a_min)/n

p_step = 0.0015

a=a_min;
% for  p1=1.1:0.1:1.3;
while a<=a_max
    lp=FO_Lyapunov_a(4,ext_fcn,0,0.2,1,[0 0 1 1],0.2,1,a)
    a=a+p_step;
    plot(a,lp,'.')
    drawnow();
    % end
end

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

Error in solution>LE_RF_a (line 24)
X= [x(5) x(9) x(13) x(17);

Error in FDE12>f_vectorfield (line 300)
f = feval(Probl.fdefun,t,y,Probl.param) ;

Error in FDE12 (line 114)
f_temp = f_vectorfield(t0,y0(:,1),Probl) ;

Error in solution>FO_Lyapunov_a (line 61)
[T,Y] = FDE12(q,ext_fcn,t,t+h_norm,[0 0 1 1],h,a);

function f=LE_RF_a(t,x,a)
%p is the parameter
% x=1:20;
f=zeros(16,1);
b=0.4;c=11;d=6;e=130;f1=10;
X= [x(5) x(9) x(13)  x(17);
    x(6) x(10) x(14) x(18);
    x(7) x(11) x(15) x(19);
    x(8) x(12)  x(16) x(20) ];
f(1)=a*x(2)+(0.2+0.2*(abs(x(4)))).*x(3);
f(2)=b*((x(4).^2)-13).*x(3)-c*x(2);
f(3)=-d*x(1)-e*x(2)-f1*x(3);
f(4)=(x(3).^2)-x(4).^2;
Jac=[0 a 0.2+0.2*abs(x(4)) 0.2*x(4).*x(3)/abs(x(3));
    0 -c b*((x(4).^2)-13) b*((x(4)*2)-13).*x(3);
    -d -e -f1 0;
    0 0 2*x(3) -2*x(4)];

f(5:20)=Jac*X;
end
function LE=FO_Lyapunov_a(ne,ext_fcn,t_start,h_norm,t_end,x_start,h,q,a);
% Memory allocation
ne=4;
x=zeros(ne*(ne+1),1); 
x0=x;
q=1;
% y0=x_start
c=zeros(ne,1);
gsc=c; zn=c;
n_it = round((t_end-t_start)/h_norm);
% Initial values
x(1:ne)=x_start;
i=1;
while i<=ne
    x((ne+1)*i)=1.0;
    i=i+1;
end
t=t_start;
% Main loop
it=1;
while it<=n_it
    % Solutuion of extended ODE system of FO using FDE12 routine
    [T,Y] = FDE12(q,ext_fcn,t,t+h_norm,[0 0 1 1],h,a); 
    t=t+h_norm;
    Y=transpose(Y);
    x=Y(size(Y,1),:); %solution at t+h_norm
    i=1;
    while i<=ne
        j=1;
        while j<=ne;
            x0(ne*i+j)=x(ne*j+i);
            j=j+1;
        end;
        i=i+1;
    end;
    %       construct new orthonormal basis by gram-schmidt
    zn(1)=0.0;
    j=1;
    while j<=ne
        zn(1)=zn(1)+x0(ne*j+1)^2; 
        j=j+1;
    end;
    zn(1)=sqrt(zn(1));
    j=1;
    while j<=ne
        x0(ne*j+1)=x0(ne*j+1)/zn(1); 
        j=j+1;
    end
    j=2;
    while j<=ne
        k=1;
        while k<=j-1
            gsc(k)=0.0;
            l=1;
            while l<=ne;
                gsc(k)=gsc(k)+x0(ne*l+j)*x0(ne*l+k);
                l=l+1;
            end
            k=k+1;
        end
        k=1;
        while k<=ne
            l=1;
            while l<=j-1
                x0(ne*k+j)=x0(ne*k+j)-gsc(l)*x0(ne*k+l);
                l=l+1;
            end
            k=k+1;
        end;
        zn(j)=0.0;
        k=1;
        while k<=ne
            zn(j)=zn(j)+x0(ne*k+j)^2;
            k=k+1;
        end
        zn(j)=sqrt(zn(j));
        k=1;
        while k<=ne
            x0(ne*k+j)=x0(ne*k+j)/zn(j);
            k=k+1;
        end
        j=j+1;
    end
    %       update running vector magnitudes
    k=1;
    while k<=ne;
        c(k)=c(k)+log(zn(k));
        k=k+1;
    end;
    %       normalize exponent        
    k=1;
    while k<=ne
        LE(k)=c(k)/(t-t_start);
        k=k+1;
    end
    i=1;
    while i<=ne
        j=1;
        while j<=ne;
            x(ne*j+i)=x0(ne*i+j);
            j=j+1;
        end
        i=i+1;
    end;
    x=transpose(x);
    it=it+1;
end
end

[nune pratyusha]

function run_LE_FO_p1(ne,ext_fcn,t_start,h_norm,t_end,x_start,h,q,p_min,p_max,n);

figure();

hold on;

ne=3;

ext_fcn=@LE_RF_p1;

t_start=0;h_norm=0.02;t_end=1;

x_start=[0.1 0.1 0.1];h=0.2;q=1;

p_max=1.3;

p_min=1.1;

n=800;

p_step=(p_max-p_min)/n

p1=p_min;

% for p1=1.1:0.1:1.3;

while p1<=p_max

lp=FO_Lyapunov_p1(ne,ext_fcn,0,0.2,1,[0.1 0.1 0.1],0.2,1,p1)

p1=p1+p_step;

plot(p1,lp,'.')

drawnow();

% end

end

function LE=FO_Lyapunov_p1(ne,ext_fcn,t_start,h_norm,t_end,x_start,h,q,p1);

% Memory allocation

x=zeros(ne*(ne+1),1);

x0=x;

c=zeros(ne,1);

gsc=c; zn=c;

n_it = round((t_end-t_start)/h_norm);

% Initial values

x(1:ne)=x_start;

i=1;

while i<=ne

x((ne+1)*i)=1.0;

i=i+1;

end

t=t_start;

% Main loop

it=1;

while it<=n_it

% Solutuion of extended ODE system of FO using FDE12 routine

[T,Y] = FDE12(q,ext_fcn,t,t+h_norm,x,h,p1);

t=t+h_norm;

Y=transpose(Y);

x=Y(size(Y,1),:); %solution at t+h_norm

i=1;

while i<=ne

j=1;

while j<=ne;

x0(ne*i+j)=x(ne*j+i);

j=j+1;

end;

i=i+1;

end;

% construct new orthonormal basis by gram-schmidt

zn(1)=0.0;

j=1;

while j<=ne

zn(1)=zn(1)+x0(ne*j+1)^2;

j=j+1;

end;

zn(1)=sqrt(zn(1));

j=1;

while j<=ne

x0(ne*j+1)=x0(ne*j+1)/zn(1);

j=j+1;

end

j=2;

while j<=ne

k=1;

while k<=j-1

gsc(k)=0.0;

l=1;

while l<=ne;

gsc(k)=gsc(k)+x0(ne*l+j)*x0(ne*l+k);

l=l+1;

end

k=k+1;

end

k=1;

while k<=ne

l=1;

while l<=j-1

x0(ne*k+j)=x0(ne*k+j)-gsc(l)*x0(ne*k+l);

l=l+1;

end

k=k+1;

end;

zn(j)=0.0;

k=1;

while k<=ne

zn(j)=zn(j)+x0(ne*k+j)^2;

k=k+1;

end

zn(j)=sqrt(zn(j));

k=1;

while k<=ne

x0(ne*k+j)=x0(ne*k+j)/zn(j);

k=k+1;

end

j=j+1;

end

% update running vector magnitudes

k=1;

while k<=ne;

c(k)=c(k)+log(zn(k));

k=k+1;

end;

% normalize exponent

k=1;

while k<=ne

LE(k)=c(k)/(t-t_start);

k=k+1;

end

i=1;

while i<=ne

j=1;

while j<=ne;

x(ne*j+i)=x0(ne*i+j);

j=j+1;

end

i=i+1;

end;

x=transpose(x);

it=it+1;

end

function f=LE_RF_p1(t,x,p1)

%p is the parameter

f=zeros(9,1);

X= [x(4) x(7) x(10);

x(5) x(8) x(11);

x(6) x(9) x(12)];

%RF equations

f(1)=x(2).*(x(3)-1+x(1).*x(1))+0.1*x(1);

f(2)=x(1).*(3*x(3)+1-x(1).*x(1))+0.1*x(2);

f(3)=-2.*x(3).*(p1+x(1).*x(2));

%Jacobian matrix

J=[2*x(1).*x(2)+0.1, x(1).*x(1)+x(3)-1, x(2);

-3*x(1).*x(1)+3*x(3)+1,0.1,3*x(1);

-2*x(2).*x(3),-2*x(1).*x(3),-2*(x(1).*x(2)+p1)];

f(4:12)=J*X; % To be modified if ne>3

In this code i put above 4 dimensional equations and i am getting error

[nune pratyusha]

After substition 4 dimensional equations why it showing error i didn't understand

[Cris LaPierre]

What is the error? Please share the complete error message (all the red text)

[nune pratyusha]

last i shared code is running that code have no error...(3 dimensinoal)

I substitute the 4 dimensional equations it showing error like how to mention in red color i don't know that's why i mentioning error

error:

>>

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

Error in solution>LE_RF_a (line 24)

X= [x(5) x(9) x(13) x(17);

Error in FDE12>f_vectorfield (line 300)

f = feval(Probl.fdefun,t,y,Probl.param) ;

Error in FDE12 (line 114)

f_temp = f_vectorfield(t0,y0(:,1),Probl) ;

Error in solution>FO_Lyapunov_a (line 61)

[T,Y] = FDE12(q,ext_fcn,t,t+h_norm,[0 0 1 1],h,a);

[nune pratyusha]

error

[Cris LaPierre]

I have already explained this error in a previous comment.

To summarize, in the code you shared for your 3D equation, you create the x input dynamically and pass that into FDE12.

ne=3;

...

x=zeros(ne*(ne+1),1); % 12x1

...

[T,Y] = FDE12(q,ext_fcn,t,t+h_norm,x,h,p1);

in the code you shared for your 4D equation, you also dynamically create x, but then hardcode your input as a 1x4 vector.

ne=4;

x=zeros(ne*(ne+1),1); % 20x1

...

[T,Y] = FDE12(q,ext_fcn,t,t+h_norm,[0 0 1 1],h,a);

Use the same syntax in both cases to obtain the same behavior

[nune pratyusha]

it's working thank you so much