Runge Kutta 4 method to solve second order ODE

Please help. I have been stuck at it for a while:

I am trying to solve a second order differential equation where U_dot= V and V_dot = d*U-c*U^3-b*V+a*sin(w*t)

Now I made the code (analytical part is meant for double checking myself, ignore it). My code gives an error saying "index exceeds array bounds". The loop refuses to accept anything other than for i = 1:1. How can I run this loop with Runge Kutta calculations?

Thank you!

w = 1.3;
dt = 2*pi/(w*100);
a= 0.25;
b=0.1;
c=1;
d = 1;
x0=1;
y0=0;
t=0:dt:5000;
transient = 4250;
tran_str= 300;
% %analytical
% w0=sqrt(-d);
% A=sqrt(x0^2+(y0/w0)^2);
% phi = atan(x0*w0/y0);
% xan = A*sin(w0*t+phi);
% yan = A*w0*cos(w0*t+phi);
%RK4
f1 = @(t,x,y) y;
f2 = @(t,x,y) d*x-c*x^3-b*y+a*sin(w*t);
h=dt
for i=1:length(t-1)
  t(i+1) = t(i)+i*h;
  k1x = f1(t(i),x0(i),y0(i));
  k2x = f2(t(i)+0.5*h,x0(i)+0.5*k1y*h, y0(i)+0.5*k1y*h);
  k3x = f2(t(i)+0.5*h,x0(i)+0.5*k2y*h, y0(i)+0.5*k2y*h);
  k4x = f2(t(i)+0.5*h,x0(i)+0.5*k3y*h, y0(i)+0.5*k3y*h);
  
  k1y = f2(t(i),x0(i),y0(i));
  k2y = f2(t(i)+0.5*h,x0(i)+0.5*k1y*h, y0(i)+0.5*k1y*h);
  k3y = f2(t(i)+0.5*h,x0(i)+0.5*k2y*h, y0(i)+0.5*k2y*h);
  k4y = f2(t(i)+0.5*h,x0(i)+0.5*k3y*h, y0(i)+0.5*k3y*h);
  y(i+1) = y0(i)+1/6*(k1y+2*k2y+2*k3y+k4y);  
    
    


end 

[Haya M]

Do you use your f1 inside the loop?

[Alina Li]

Thank you for the question, it actually made me realize that I should have included it. I edited the code above. Unfrotunately, the same error persists.

[Haya M]

What do you mean by k1y in k2x and check the rest as well.

Answer 1

[Alan Stevens]

Your integration loop is mightily scrambled. It should be like this

x(1) = x0;
y(1) = y0;
for i=1:length(t)-1
  t(i+1) = i*h;  
  k1x = f1(t(i),x(i),y(i));
  k1y = f2(t(i),x(i),y(i));
  
  k2x = f1(t(i)+0.5*h,x(i)+0.5*k1x*h, y(i)+0.5*k1y*h);
  k2y = f2(t(i)+0.5*h,x(i)+0.5*k1x*h, y(i)+0.5*k1y*h);
  
  k3x = f1(t(i)+0.5*h,x(i)+0.5*k2x*h, y(i)+0.5*k2y*h);
  k3y = f2(t(i)+0.5*h,x(i)+0.5*k2x*h, y(i)+0.5*k2y*h);
  
  k4x = f1(t(i)+h,x(i)+k3x*h, y(i)+k3y*h);
  k4y = f2(t(i)+h,x(i)+k3x*h, y(i)+k3y*h);
  
  x(i+1) = x(i)+1/6*(k1x+2*k2x+2*k3x+k4x)*h; 
  y(i+1) = y(i)+1/6*(k1y+2*k2y+2*k3y+k4y)*h;  
    
end 

(However, I don't think your analytical solution is the solution to your equations.)

[Bruce Taylor]

This is a beautiful solution. I modified it to analyze a series R-L-C circuit and compared the result to the exact solution...perfect match.

Bruce Taylor