RK4 for a system of 2 second order DEs

[aoeuDvorakhtn]

1. For this problem, we are asked to design a MatLab/Octave code that will calculate and plot the orbit of a sattelite



2. We are given the following 2nd order equations
x''= \frac{-x}{( \sqrt{x^{2}+y^{2}})^{3}}
y''= \frac{-y}{( \sqrt{x^{2}+y^{2}})^{3}}

We can convert this into the following system of 4 1st order DEs by substitiuting z=x' and w=y'
x'=z
y'=w
z'= \frac{-x}{( \sqrt{x^{2}+y^{2}})^{3}}
w'= \frac{-y}{( \sqrt{x^{2}+y^{2}})^{3}}

We use the following initial vallues
[x0 y0 z0 w0] = [4 0 0 0.5]

The Attempt at a Solution

I managed to build a program that works well enough that, after adding annotations, I could hand it in. However, as you can see below, the lines within the for loop could use a little help as far as elegance goes (I've attached the .m file if my copy/paste is too difficult to read). Ideally, I'd like to be able to make those lines a little more compact. Any ideas?

X0=[4 0 0 0.5];
t0=0;
tmax=50;
h=0.25;
SOL=[t0 X0;];
f1=inline('z','z');
f2=inline('w','w');
f3=inline('-x/((sqrt(x^2+y^2))^3)','x','y');
f4=inline('-y/((sqrt(x^2+y^2))^3)','x','y');
for i = 1:tmax/h
g1=[f1(X0(3)) f2(X0(4)) f3(X0(1),X0(2)) f4(X0(1),X0(2))];
g2=[f1(X0(3)+0.5*h*g1(3)) f2(X0(4)+0.5*h*g1(4)) f3(X0(1)+0.5*h*g1(1),X0(2)+0.5*h*g1(2)) f4(X0(1)+0.5*h*g1(1),X0(2)+0.5*h*g1(2))];
g3=[f1(X0(3)+0.5*h*g2(3)) f2(X0(4)+0.5*h*g2(4)) f3(X0(1)+0.5*h*g2(1),X0(2)+0.5*h*g2(2)) f4(X0(1)+0.5*h*g2(1),X0(2)+0.5*h*g2(2))];
g4=[f1(X0(3)+h*g3(3)) f2(X0(4)+h*g3(4)) f3(X0(1)+h*g3(1),X0(2)+h*g3(2)) f4(X0(1)+h*g3(1),X0(2)+h*g3(2))];
X1=X0+(h/6)*(g1+2*g2+2*g3+g4);
t1=t0+h;
SOL=[SOL;t1 X1];
X0=X1;
t0=t1;
end
plot(SOL(:,2),SOL(:,3))
axis('equal')

 

[aoeuDvorakhtn]

In case anyone is interested, I did figure something out. It added a few lines to the code, but it did make the loop a bit prettier.

%sets initial position and velocity values for x and y
%(z and w are the velecities of x and y, respectively)
X0=[4 0 0 0.5];
%initializes the time variable (x, y, z, and w are functions of time)
%and sets a maximum t for which to approximate
t0=0;
tmax=50;
%stepsize
h=0.25;
%initializes array in which to store approximated valies.
SOL=[t0 X0];
%these inline functions are our system of 4 1st order DEs, x', y', z', and w'
f1=inline('z','z');
f2=inline('w','w');
f3=inline('-x/((x^2+y^2)^1.5)','x','y');
f4=inline('-y/((x^2+y^2)^1.5)','x','y');
%approximates x, y, z, and w in increments of 0.0025 from t=0 to t=50 using the
%classical 4th order Runge-Kutta method
for i = 1:tmax/h

g1=[f1(X0(3)) f2(X0(4)) f3(X0(1),X0(2)) f4(X0(1),X0(2))];
k1=X0+0.5*h*g1;
g2=[f1(k1(3)) f2(k1(4)) f3(k1(1),k1(2)) f4(k1(1),k1(2))];
k2=X0+0.5*h*g2;
g3=[f1(k2(3)) f2(k2(4)) f3(k2(1),k2(2)) f4(k2(1),k2(2))];
k3=X0+h*g3;
g4=[f1(k3(3)) f2(k3(4)) f3(k3(1),k3(2)) f4(k3(1),k3(2))];
X1=X0+(h/6)*(g1+2*g2+2*g3+g4); %RK4
t1=t0+h;
SOL=[SOL;t1 X1];
X0=X1;
t0=t1;​

end
%plots our approximation to the position functions, x and y and sets the axes
%to use equal scales
plot(SOL(:,2),SOL(:,3))
axis('equal')