MATLAB: Runge Kutta method for coupled oscillator system.

coupled oscillator

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));

Best Answer

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')
Related Question