Index exceeds the number of array elements (1)

Can anyone please help me to resolve my issue. I have been trying to resolve it since two days. It's not easy for me as I am new to MATLAB.

close all; clear all; clc;
%initializing x,y,t
h=0.01;   %step size
t=0:h:300;
x1 = zeros(1,numel(t)); y1 = x1; x2 = x1; y2 = x1;
x1(1)=1;
y1(1)=1;
x2(1)=2;
y2(1)=2;
%value of constants
a=0.1;
b=0.3;
omega=4;
Cnp=0.2;
G=1;
%ode
p=@(x1,y1,x2,y2) (a-x1^2-y1^2)*x1-omega*y1+G*Cnp*(x2-x1);
q=@(x1,y1,x2,y2) (a-x1^2-y1^2)*y1+omega*x1+G*Cnp*(y2-y1);
%loop
for i=1:(length(t)-1)
    k1=p(x1(i),y1(i),x2(i),y2(i));
    l1=q(x1(i),y1(i),x2(i),y2(i));
 
    
    k2=p((x1(i)+(h/2)*k1),(y1(i)+(h/2)*l1),(x2(i)+(h/2)*k1),(y2(i)+(h/2)*l1));
    l2=q((x1(i)+(h/2)*k1),(y1(i)+(h/2)*l1),(x2(i)+(h/2)*k1),(y2(i)+(h/2)*l1));
      
    k3=p((x1(i)+(h/2)*k2),(y1(i)+(h/2)*l2),(x2(i)+(h/2)*k2),(y2(i)+(h/2)*l2));
    l3=q((x1(i)+(h/2)*k2),(y1(i)+(h/2)*l2),(x2(i)+(h/2)*k2),(y2(i)+(h/2)*l2));
      
    k4=p((x1(i)+k3*h),(y1(i)+l3*h),(x2(i)+k3*h),(y2(i)+l3*h));
    l4=q((x1(i)+k3*h),(y1(i)+l3*h),(x2(i)+k3*h),(y2(i)+l3*h));
    
      x1(i+1) = x1(i) + h*(k1+2*k2+2*k3+k4)/6;
      y1(i+1) = y1(i) + h*(l1+2*l2+2*l3+l4)/6;
      x2(i+1) = x2(i) + h*(k1+2*k2+2*k3+k4)/6;
      y2(i+1) = y2(i) + h*(l1+2*l2+2*l3+l4)/6;
end
%draw
plot(t,x1,'r',t,x2,'y')
xlabel('t','fontsize',14,'fontweight','bold')
ylabel('x1 & x2','fontsize',14,'fontweight','bold')
set(gca,'Color','k')
legend('x1','x2','TextColor','w')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
x1=x1(101:300);
y1=y1(101:300);
x2=x2(101:300);
y2=y2(101:300);
X1dot=(x1(2:101)-x1(1:100))/h;
X2dot=(x2(2:101)-x2(1:100))/h;
Y1dot=(y1(2:101)-y1(1:100))/h;
Y2dot=(y2(2:101)-y2(1:100))/h;
hold on

M=zeros([100,10]);
for i=1:100
     for j=1:32
         if j==1
             M(i,j)=1;
            elseif j==2
               M(i,j)=x1(i);
            elseif j==3
               M(i,j)=x2((i));
            elseif j==4
               M(i,j)=y1((i));
            elseif j==5
               M(i,j)=y2((i));
            elseif j==6
               M(i,j)=x1(i)*x2(i);
            elseif j==7
                M(i,j)=y1((i))*y2((i));
            elseif j==8
                M(i,j)=x1((i))*y1((i));
            elseif j==9
                M(i,j)=x1((i))*y2((i));
            elseif j==10
                M(i,j)=x2((i))*y1((i));
            elseif j==11
                M(i,j)=x2((i))*y2((i));   
            elseif j==12
                M(i,j)=x1(i)*x2(i)*y1(i);
            elseif j==13
                M(i,j)=x1(i)*x2(i)*y2(i);
            elseif j==14
                M(i,j)=x1(i)*y1(i)*y2(i); 
            elseif j==15
                M(i,j)=x2(i)*y1(i)*y2(i);  
            elseif j==16
                M(i,j)=x1(i)*x2(i)*y1(i)*y2(i);  
            elseif j==17
                M(i,j)=x1(i)^2;
            elseif j==18
                M(i,j)=x2(i)^2;
            elseif j==19
                M(i,j)=y1(i)^2;
            elseif j==20
                M(i,j)=y2(i)^2;
            elseif j==21
                M(i,j)=x1(i)^2*x2(i);
            elseif j==22
                M(i,j)=x1(i)^2*y1(i);
            elseif j==23
                M(i,j)=x1(i)^2*y2(i);
            elseif j==24
                M(i,j)=x2(i)^2*x1(i);
            elseif j==25
                M(i,j)=x2(i)^2*y1(i);     
            elseif j==26
                M(i,j)=x2(i)^2*y2(i);
            elseif j==27
                M(i,j)=y1(i)^2*x1(i);
            elseif j==28
                M(i,j)=y1(i)^2*x2(i);
            elseif j==29
                M(i,j)=y1(i)^2*y2(i);
            elseif j==30
                M(i,j)=y2(i)^2*x1(i);
            elseif j==31
                M(i,j)=y2(i)^2*x2(i);
            else
                 M(i,j)=y2(i)^2*y1(i);     
            end
      end
end      
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                 Reconstructed Signal                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x1_rec(1)=x1(1);
y1_rec(1)=y1(1);
x2_rec(1)=x2(1);
y2_rec(1)=y2(1);
%value of constants
%a1=x1
a1=0.0305;
a2=0.0000;
a3=0.0000;
a4=-3.2154;
a5=0.0000;
a6=0.0000;
a7=0.0000;
a8=0.0000;
a9=0.0000;
a10=-0.5857;
a11=0.0000;
a12=0.0000;
a13=0.0000;
a14=0.0000;
a15=0.0000;
a16=0.0000;
a17=0.0000;
a18=0.0000;
a19=0.0000;
a20= 0.2412;
a21=0.0000;
a22=0.0000;
a23=0.0000;
a24=0.0000;
a25=0.0000;
a26=0.0000;
a27=0.0000;
a28=0.0000;
a29=0.0000;
a30=0.0000;
a31=0.0000;
a32=0.0000;
%c=x2
c1=0.0305;
c2=0.0000;
c3=0.0000;
c4=-3.2154;
c5=0.0000;
c6=0.0000;
c7=0.0000;
c8=0.0000;
c9=0.0000;
c10=-0.5857;
c11=0.0000;
c12=0.0000;
c13=0.0000;
c14=0.0000;
c15=0.0000;
c16=0.0000;
c17=0.0000;
c18=0.0000;
c19=0.0000;
c20=0.2412;
c21=0.0000;
c22=0.0000;
c23=0.0000;
c24=0.0000;
c25=0.0000;
c26=0.0000;
c27=0.0000;
c28=0.0000;
c29=0.0000;
c30=0.0000;
c31=0.0000;
c32=0.0000;
h=0.01;   %step size
t=0:h:200;
%ode
p=@(x1_rec,y1_rec,x2_rec,y2_rec) a1+a2*x1_rec+a3*x2_rec+a4*y1_rec+a5*y2_rec+a6*x1_rec*x2_rec+a7*y1_rec*y2_rec+a8*x1_rec*y1_rec+a9*x1_rec*y2_rec+a10*x2_rec*y1_rec+a11*x2_rec*y2_rec+a12*x1_rec*x2_rec*y1_rec+a13*x1_rec*x2_rec*y2_rec+a14*x1_rec*y1_rec*y2_rec+a15*x2_rec*y1_rec*y2_rec+a16*x1_rec*x2_rec*y1_rec*y2_rec+a17*x1_rec^2+a18*x2_rec^2+a19*y1_rec^2+a20*y2_rec^2+a21*x1_rec^2*x2_rec+a22*x1_rec^2*y1_rec+a23*x1_rec^2*y2_rec+a24*x2_rec^2*x1_rec+a25*x2_rec^2*y1_rec+a26*x2_rec^2*y2_rec+a27*y1_rec^2*x1_rec+a28*y1_rec^2*x2_rec+a29*y1_rec^2*y2_rec+a30*y2_rec^2*x1_rec+a31*y2_rec^2*x2_rec+a32*y2_rec^2*y1_rec;
q=@(x1_rec,y1_rec,x2_rec,y2_rec) c1+c2*x1_rec+c3*x2_rec+c4*y1_rec+c5*y2_rec+c6*x1_rec*x2_rec+c7*y1_rec*y2_rec+c8*x1_rec*y1_rec+c9*x1_rec*y2_rec+c10*x2_rec*y1_rec+c11*x2_rec*y2_rec+c12*x1_rec*x2_rec*y1_rec+c13*x1_rec*x2_rec*y2_rec+c14*x1_rec*y1_rec*y2_rec+c15*x2_rec*y1_rec*y2_rec+c16*x1_rec*x2_rec*y1_rec*y2_rec+c17*x1_rec^2+c18*x2_rec^2+c19*y1_rec^2+c20*y2_rec^2+c21*x1_rec^2*x2_rec+c22*x1_rec^2*y1_rec+c23*x1_rec^2*y2_rec+c24*x2_rec^2*x1_rec+c25*x2_rec^2*y1_rec+c26*x2_rec^2*y2_rec+c27*y1_rec^2*x1_rec+c28*y1_rec^2*x2_rec+c29*y1_rec^2*y2_rec+c30*y2_rec^2*x1_rec+c31*y2_rec^2*x2_rec+c32*y2_rec^2*y1_rec;
%loop
for i=1:(length(t)-1)
    k1=p(x1_rec(i),y1_rec(i),x2_rec(i),y2_rec(i));
    l1=q(x1_rec(i),y1_rec(i),x2_rec(i),y2_rec(i));
 
    
    k2=p((x1_rec(i)+(h/2)*k1),(y1_rec(i)+(h/2)*l1),(x2_rec(i)+(h/2)*k1),(y2_rec(i)+(h/2)*l1));
    l2=q((x1_rec(i)+(h/2)*k1),(y1_rec(i)+(h/2)*l1),(x2_rec(i)+(h/2)*k1),(y2_rec(i)+(h/2)*l1));
      
    k3=p((x1_rec(i)+(h/2)*k2),(y1_rec(i)+(h/2)*l2),(x2_rec(i)+(h/2)*k2),(y2_rec(i)+(h/2)*l2));
    l3=q((x1_rec(i)+(h/2)*k2),(y1_rec(i)+(h/2)*l2),(x2_rec(i)+(h/2)*k2),(y2_rec(i)+(h/2)*l2));
      
    k4=p((x1_rec(i)+k3*h),(y1_rec(i)+l3*h),(x2_rec(i)+k3*h),(y2_rec(i)+l3*h));
    l4=q((x1_rec(i)+k3*h),(y1_rec(i)+l3*h),(x2_rec(i)+k3*h),(y2_rec(i)+l3*h));
    
      x1_rec(i+1) = x1_rec(i) + h*(k1+2*k2+2*k3+k4)/6;
      x2_rec(i+1) = x2_rec(i) + h*(k1+2*k2+2*k3+k4)/6;
end

plot(x1_rec(1:length(x1)),'r','LineWidth',2)
plot(x2_rec(1:length(x2)),'b','LineWidth',2)



%%%%%Error%%%%

Index exceeds the number of array elements (1).

Error in coup_withoutcvx (line 208)
    k1=p(x1_rec(i),y1_rec(i),x2_rec(i),y2_rec(i));

Answer 1

[Walter Roberson]

You increase the size of x1_rec and x2_rec in the loop, but you do not increase the size of y2_rec, so when i becomes 2 then

    k1=p(x1_rec(i),y1_rec(i),x2_rec(i),y2_rec(i));

is going to require y2_rec(2) which does not exist.

[Heya :)]

Where do I have to change the code to increase the size of y2_rec?

[Walter Roberson]

Right after

      x2_rec(i+1) = x2_rec(i) + h*(k1+2*k2+2*k3+k4)/6;

assign to y1_rec(i+1) and y2_rec(i+1)

However I do not know what values should be assigned to them.

[Heya :)]

Can you alter the code please?

[Walter Roberson]

x2_rec(i+1) = x2_rec(i) + h*(k1+2*k2+2*k3+k4)/6;

y1_rec(i+1) = rand;
y2_rec(i+1) = rand;

... since I have no idea what the formulas should be.

Your code looks like it might be a Runge-Kutta method. If so you should probably study the material on Runge-Kutta.

Answer 2

[Geoff Hayes]

The error is here

    k1=p(x1_rec(i),y1_rec(i),x2_rec(i),y2_rec(i));
    l1=q(x1_rec(i),y1_rec(i),x2_rec(i),y2_rec(i));

where the y1_rec and the y2_rec are 1x1 scalars yet you are treating them as if they are arrays. Since i is the for loop "step" variable, then on the second iteration of this loop you will get the "index exceeds number of array elements". You either need to treat them as scalars or update the y1_rec and y2_rec arrays on each iteration of the loop like you do for the x1_rec and x2_rec.

[Heya :)]

How to make y1_rec and y2_rec arrays?