MATLAB: RK4 orbit problem

arraysorbitsrk4

% Parameters
h = 0.0027;      %step size
tfinal = 5;  % t=[0,tfinal]
GM = 4*pi^2 ;
% Initial Conditions
vx = zeros(1,1825);
vy = zeros(1,1825);
x = zeros(1,1825);
y = zeros(1,1825);
x(1) = 0;
y(1) = 1;
vy(1) = 0;
vx(1) = 1;
t = zeros(1,1825);
% Define the ODE function handle
dvxdt=@(t,x,vx) GM*x./(x.^2+y.^2).^3;
dvydt=@(t,y,vy) GM*y./(x.^2+y.^2).^3;
dxdt=@(t,x,vx) vx;
dydt=@(t,y,vy) vy;
%Runge Kutta Loop
for i=1:ceil(tfinal/h)
    t(i+1) = t(i) + h;
    
    
    k1vx = dvxdt(t(i)     ,x(i)     ,vx(i));
    k1vy = dvydt(t(i)     ,y(i)     ,vy(i));
    k1x = dxdt(t(i)     ,x(i)     ,vx(i));
    k1y = dydt(t(i)     ,y(i)     ,vy(i));
    
    k2vx = dvxdt(t(i)+0.5*h,x(i)+0.5*h*k1vx,vx(i)+0.5*h*k1vx);
    k2vy = dvydt(t(i)+0.5*h,y(i)+0.5*h*k1vy,vy(i)+0.5*h*k1vy);
    k2x = dxdt(t(i)+0.5*h,x(i)+0.5*h*k1x,vx(i)+0.5*h*k1x);
    k2y = dydt(t(i)+0.5*h,y(i)+0.5*h*k1y,vy(i)+0.5*h*k1y);
    
    k3vx = dvxdt(t(i)+0.5*h,x(i)+0.5*h*k2vx,vx(i)+0.5*h*k2vx);
    k3vy = dvydt(t(i)+0.5*h,y(i)+0.5*h*k2vy,vy(i)+0.5*h*k2vy);
    k3x = dxdt(t(i)+0.5*h,x(i)+0.5*h*k2x,vx(i)+0.5*h*k2x);
    k3y = dydt(t(i)+0.5*h,y(i)+0.5*h*k2y,vy(i)+0.5*h*k2y);
    
    k4vx = dvxdt(t(i)+    h,x(i)    +h*k3vx,vx(i)     +h*k3vx);
    k4vy = dvydt(t(i)+    h,y(i)    +h*k3vy,vy(i)     +h*k3vy);
    k4x = dxdt(t(i)+    h,x(i)    +h*k3x,vx(i)     +h*k3x);
    k4y = dydt(t(i)+    h,y(i)    +h*k3y,vy(i)     +h*k3y);
    
    vx(i+1)=vx(i)+h/6*(k1vx+2*k2vx+2*k3vx+k4vx);
    vy(i+1)=vy(i)+h/6*(k1vy+2*k2vy+2*k3vy+k4vy);
    x(i+1)=x(i)+h/6*(k1x+2*k2x+2*k3x+k4x);
    y(i+1)=y(i)+h/6*(k1y+2*k2y+2*k3y+k4y);
end
%plot the results
plot(x,y)
it says the left and right sides have a different number of elements and I dont know how to fix it.

Best Answer

dvxdt=@(t,x,vx) GM*x./(x.^2+y.^2).^3;
that accepts vx as a parameter but does not use it.
It does not have any parameter named y but uses y. In that case, it is going to use the variable value that the variable y had at the time the anonymous function was created, which comes from
y = zeros(1,1825);
This is a vector, and even though it only contributes 0 at all components, it makes the result of dvxdt be a vector. That becomes a problem because it means that
    k1vx = dvxdt(t(i)     ,x(i)     ,vx(i));
gives a vector result, so
vx(i+1)=vx(i)+h/6*(k1vx+2*k2vx+2*k3vx+k4vx);
has a vector on the right hand side, but designates a scalar output location.
Related Question