Error in plotting the waves.

I am supposed to get the two waves as showing in the picture with this code. Can anybody please spot the error or refine the code

as I am getting two straight lines only.

close all; clear all; clc;
%initializing x,y,t
h=0.05;   %step size
t=0:h:400;
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                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
h=0.05;   %step size
t=0:h:400;
x1_rec = zeros(1,numel(t));
y1_rec = x1_rec;
x2_rec = x1_rec;
y2_rec = x1_rec;
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;
%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=2:test_points
    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;
      %y1(i+1) = y1(i) + h*(l1+2*l2+2*l3+l4)/6;
      x2_rec(i+1) = x2_rec(i) + h*(k1+2*k2+2*k3+k4)/6;
      %y2(i+1) = y2(i) + h*(l1+2*l2+2*l3+l4)/6;
end

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

Answer 1

[Walter Roberson]

You initialize y1_rec and y2_rec to all zero, and you assign 0 to y1_rec(1) and y2_rec(1), but you never assign anything non-zero to the entries.

[Heya :)]

How to correct it? Please share any command etc or a code.

[Heya :)]

If I don't use this command [x1_rec = zeros(1,numel(t));] then the size don't match.

[Walter Roberson]

We went through this already in https://www.mathworks.com/matlabcentral/answers/662703-index-exceeds-the-number-of-array-elements-1#comment_1161423

To summarize: you did not document the purpose of the code, and you did not provide a link to a description of the algorithm, so none of us has any idea what y1_rec and y2_rec should be.