Manual Runge-Kutta 2 for a system of 4 ODE's

Question

A follow-up on this question :

I have the following system:

$\frac{dx}{dt} = p_x \\ \frac{dy}{dt} = p_y \\ \frac{dp_x}{dt}=-\frac{\partial V}{\partial x} \\ \frac{dp_y}{dt}=-\frac{\partial V}{\partial y}$

,where $V(x,y) = \frac{1}{2} (\omega_x x^2 + \omega_y y^2)$, and $\omega_x = 1, \omega_y = 2$. Following the advice in the linked page, I wrote the system in vector form, and now I would like to manually run a Do loop in order to solve it with a Runge-Kutta 2 method. The loop itself isn't difficult, and before the loop are the initial data

tf=100;
dt=1/10;
T=Table[i,{i,ti,tf,dt}];
X=ConstantArray[0,Length[T]];
Y=ConstantArray[0,Length[T]];
PX=ConstantArray[0,Length[T]];
PY=ConstantArray[0,Length[T]];
PX[[1]]=1;
PY[[1]]=1;
Do[k1=dt*f[T[[n]],X[[n]],Y[[n]],PX[[n]],PY[[n]]];
l1=dt*g[T[[n]],X[[n]],Y[[n]],PX[[n]],PY[[n]]];
m1=dt*h[T[[n]],X[[n]],Y[[n]],PX[[n]],PY[[n]]];
n1=dt*p[T[[n]],X[[n]],Y[[n]],PX[[n]],PY[[n]]];
k2=dt*f[(T[[n]]+dt/2),(X[[n]]+k1/2),(Y[[n]]+l1/2),(PX[[n]]+m1/2),(PY[[n]]+n1/2)];
l2=dt*g[(T[[n]]+dt/2),(X[[n]]+k1/2),(Y[[n]]+l1/2),(PX[[n]]+m1/2),(PY[[n]]+n1/2)];
m2=dt*h[(T[[n]]+dt/2),(X[[n]]+k1/2),(Y[[n]]+l1/2),(PX[[n]]+m1/2),(PY[[n]]+n1/2)];
n2=dt*p[(T[[n]]+dt/2),(X[[n]]+k1/2),(Y[[n]]+l1/2),(PX[[n]]+m1/2),(PY[[n]]+n1/2)];
X[[n+1]]=X[[n]]+k2;
Y[[n+1]]=Y[[n]]+l2;
PX[[n+1]]=PX[[n]]+m2;
PX[[n+1]]=PY[[n]]+n2, {n,1, Length[X]-1}]

However, I'm having trouble understanding how to apply the code to my vector of equations: How should I define the functions $f,g,h$ and $p$?

Edit: Thanks to everyone who contributed. I also came up with a simple solution to define the function in vector form, which also shortens the code considerably, so I'm uploading it for anyone who might be interested:

n1=1000;
h=0.1;
T=Table[i,{i,0,100,h}];
F[t_,{x_,y_,px_,py_}]:={px,py,-x,-2*y};
U=ConstantArray[{0.,0.,1.,1.},n1+1];

Do[K1=h*F[T[[i]],U[[i]]];   
K2=h*F[T[[i]]+h/2,U[[i]]+K1/2];
U[[i+1]]=U[[i]]+K2,
{i,1,n1}]

X=U[[All,1]];
Y=U[[All,2]];
PX=U[[All,3]];
PY=U[[All,4]];
ListPlot[Transpose[{T,X}],PlotLabel->"x(t)",AxesLabel->{"t","x"}]
ListPlot[Transpose[{T,Y}],PlotLabel->"y(t)",AxesLabel->{"t","y"}]

energy=(PX^2 + PY^2)/2 + (X^2 + 2*(Y^2))/2;
ListPlot[Transpose[{T,energy}],PlotLabel->"Energy",AxesLabel->{"t","E(t)"}]
ListPlot[Transpose[{T,energy-energy[[1]]}],PlotLabel->"Energy error",AxesLabel->{"t","E(t)-E(0)"}]

Answer 1

Using JBuck's code with implementation the midpoint method of rk2 family we have

Clear["Global`*"]

V = 1/2 (x1^2 + 2 y1^2);

ti = 0; tf = 10;
dt = 1./10;
T = Table[i, {i, ti, tf, dt}];
X = ConstantArray[0, Length[T]];
Y = ConstantArray[0, Length[T]];
PX = ConstantArray[0, Length[T]];
PY = ConstantArray[0, Length[T]];
PX[[1]] = 1;
PY[[1]] = 1;
f[t_, x_, y_, px_, py_] := px; g[t_, x_, y_, px_, py_] := py; 
h[t_, x_, y_, px_, py_] := -D[V, x1] /. {x1 -> x, y1 -> y}; 
p[t_, x_, y_, px_, py_] := -D[V, y1] /. {x1 -> x, y1 -> y};
Do[k1 = dt f[T[[n]], X[[n]], Y[[n]], PX[[n]], PY[[n]]];
 l1 = dt g[T[[n]], X[[n]], Y[[n]], PX[[n]], PY[[n]]];
 m1 = dt h[T[[n]], X[[n]], Y[[n]], PX[[n]], PY[[n]]];
 n1 = dt p[T[[n]], X[[n]], Y[[n]], PX[[n]], PY[[n]]];
 k2 = dt f[T[[n]] + dt/2, X[[n]] + k1/2, Y[[n]] + l1/2, 
    PX[[n]] + m1/2, PY[[n]] + n1/2];
 l2 = dt g[T[[n]] + dt/2, X[[n]] + k1/2, Y[[n]] + l1/2, 
    PX[[n]] + m1/2, PY[[n]] + n1/2];
 m2 = dt h[T[[n]] + dt/2, X[[n]] + k1/2, Y[[n]] + l1/2, 
    PX[[n]] + m1/2, PY[[n]] + n1/2];
 n2 = dt p[T[[n]] + dt/2, X[[n]] + k1/2, Y[[n]] + l1/2, 
    PX[[n]] + m1/2, PY[[n]] + n1/2];
 X[[n + 1]] = X[[n]] + k2;
 Y[[n + 1]] = Y[[n]] + l2;
 PX[[n + 1]] = PX[[n]] + m2;
 PY[[n + 1]] = PY[[n]] + n2;, {n, 1, Length[T] - 1}]

Visualization of numerical solution (red points) with NDSolve solution (solid lines)

sol = NDSolve[{x''[t] == -D[V, x1] /. {x1 -> x[t], y1 -> y[t]}, 
   y''[t] == -D[V, y1] /. {x1 -> x[t], y1 -> y[t]}, x[0] == 0, 
   y[0] == 0, x'[0] == 1, y'[0] == 1}, {x, y}, {t, ti, tf}];

{Show[Plot[x[t] /. sol[[1]], {t, ti, tf}, PlotLabel -> "X", 
   AxesLabel -> Automatic], 
  ListPlot[Transpose[{T, X}], PlotRange -> All, PlotStyle -> Red]], 
 Show[Plot[y[t] /. sol[[1]], {t, ti, tf}, PlotLabel -> "Y", 
   AxesLabel -> Automatic], 
  ListPlot[Transpose[{T, Y}], PlotRange -> All, PlotStyle -> Red]]} 

Update 1. We can organize rk2 step in separate module as follows

Clear["Global`*"]
U = 1/2 (x1^2 + 2 x2^2); x = ConstantArray[0, 4]; 
f[t_, x_] := {x[[3]], 
  x[[4]], -D[U, x1] /. {x1 -> x[[1]], 
    x2 -> x[[2]]}, -D[U, x2] /. {x1 -> x[[1]], x2 -> x[[2]]}};

rk2[f_, h_][{t_, x_}] := Module[{k1, k2}, k1 = f[t, x];
  k2 = f[t + h/2, x + h k1/2];
  {t + h, x + h k2}]

We can use NestList instead of Do loop and compare numerical solution with NDSolve as well

tf = 100; dt = 1/10; sol = 
 NestList[rk2[f, dt], {0, {0, 0, 1, 1}}, Round[tf/dt]];

sol1 = NDSolve[{X''[T] == -D[U, x1] /. {x1 -> X[T], x2 -> Y[T]}, 
   Y''[T] == -D[U, x2] /. {x1 -> X[T], x2 -> Y[T]}, X[0] == 0, 
   Y[0] == 0, X'[0] == 1, Y'[0] == 1}, {X, Y}, {T, 0, tf}]; {Show[
  Plot[X[T] /. sol1[[1]], {T, 0, tf}, PlotLabel -> "X"], 
  ListPlot[Transpose[{sol[[All, 1]], sol[[All, 2, 1]]}], 
   PlotStyle -> Red]], 
 Show[Plot[Y[T] /. sol1[[1]], {T, 0, tf}, PlotLabel -> "Y"], 
  ListPlot[Transpose[{sol[[All, 1]], sol[[All, 2, 2]]}], 
   PlotStyle -> Red]]}