Runge Kutta method for coupled oscillator system.

I am trying to solve these equations with the help of Runge Kutta Method. I am not sure whether I am writing the code correctly or we can use this method for coupled also getting this error (mentioned at the end of the code). Please help me to refine my code.

close all; clear all; clc;
%initializing x,y,t
t(1)=0;
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;
h=0.1;   %step size
t=0:h:50;
%ode
p=@(t,x1,y1,x2,y2) (a-x1^2-y1^2)*x1-omega*y1+G*Cnp(x2-x1);
q=@(t,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(t(i),x1(i),y1(i),x2(i),y2(i));
    l1=q(t(i),x1(i),y1(i),x2(i),y2(i));
 
    
    k2=p(t(i)+h/2,(x1(i)+(h/2)*k1),(y1(i)+(h/2)*l1),(x2(i)+(h/2)*k1),(y2(i)+(h/2)*l1));
    l2=q(t(i)+h/2,(x1(i)+(h/2)*k1),(y1(i)+(h/2)*l1),(x2(i)+(h/2)*k1),(y2(i)+(h/2)*l1));
      
    k3=p(t(i)+h/2,(x1(i)+(h/2)*k2),(y1(i)+(h/2)*l2),(x2(i)+(h/2)*k2),(y2(i)+(h/2)*l2));
    l3=q(t(i)+h/2,(x1(i)+(h/2)*k2),(y1(i)+(h/2)*l2),(x2(i)+(h/2)*k2),(y2(i)+(h/2)*l2));
      
    k4=p(t(i)+h,(x1(i)+k3*h),(y1(i)+l3*h),(x2(i)+k3*h),(y2(i)+l3*h));
    l4=q(t(i)+h,(x1(i)+k3*h),(y1(i)+l3*h),(x2(i)+k3*h),(y2(i)+l3*h));
    
      x(i+1) = x(i) + h*(k1+2*k2+2*k3+k4)/6;
      y(i+1) = y(i) + h*(l1+2*l2+2*l3+l4)/6;
end
%draw
plot(t,x,'r')
xlabel('X','fontsize',14,'fontweight','bold')
ylabel('y','fontsize',14,'fontweight','bold')
set(gca,'Color','k')


****Error***


Array indices must be positive integers or logical values.

Error in coupled>@(t,x1,y1,x2,y2)(a-x1^2-y1^2)*x1-omega*y1+G*Cnp(x2-x1) (line 17)
p=@(t,x1,y1,x2,y2) (a-x1^2-y1^2)*x1-omega*y1+G*Cnp(x2-x1);

Error in coupled (line 25)
    k2=p(t(i)+h/2,(x1(i)+(h/2)*k1),(y1(i)+(h/2)*l1),(x2(i)+(h/2)*k1),(y2(i)+(h/2)*l1));

Answer 1

[Alan Stevens]

Do you mean like this

%initializing x,y,t
h=0.1;   %step size
t=0:h:50;
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')
figure
plot(t,y1,'r',t,y2,'y')
xlabel('t','fontsize',14,'fontweight','bold')
ylabel('y1 & y2','fontsize',14,'fontweight','bold')
set(gca,'Color','k')
legend('y1','y2','TextColor','w')

[Heya :)]

Thanks a lot! Can you please tell me why you use this ?

x1 = zeros(1,numel(t)); y1 = x1; x2 = x1; y2 = x1;

[Alan Stevens]

That is just setting aside storage space for those variables. It makes the loop more efficient as space doesn't then have to be found for their new elements every iteration of the loop.

[Heya :)]

Thank you :)