So, I've implemented RK4, and I'm wondering what I can do to make it more efficient? What I've got so far is below. I wish to still record all steps. I think AppendTo is doing the most damage to the time, is there a faster alternative?

rk4[f_, variables_, valtinit_, tinit_, tfinal_, nsteps_] := 
 Module[table, xlist, ylist, step, k1, k2, k3, k4, 
 xlist = tinit;
 step = N[(tfinal - tinit)/(nsteps)];
 ylist = valtinit;
 table = xlist, ylist;
 Table[
 k1 = step* f /. MapThread[Rule, variables, ylist]; (* 
 Equivalent to step* f/.Thread[Rule[variables,ylist]]*)
 k2 = step*f /. MapThread[Rule, variables, k1/2 + ylist];
 k3 = step*f /. MapThread[Rule, variables, k2/2 + ylist];
 k4 = step*f /. MapThread[Rule, variables, k3 + ylist];
 ylist += 1/6 (k1 + 2 (k2 + k3) + k4);
 xlist += step;
 AppendTo[table, xlist, ylist];
 xlist, ylist, nsteps];
 table
 ];

Example Input:

funclist = -x + y, x - y;
initials = 1, 2;
variables = x, y;
init = 0;
final = 200;
nstep = 20000;
approx = rk4[funclist, variables, initials, init, final, nstep]//AbsoluteTiming;

3.59932,...

I'd love some suggestions!

Just to give you an impression how fast things may get when you use the right tools.

For given stepsize τ and given vector field F, this creates a CompiledFunction cStep that computes a single Runge-Kutta step

F = X [Function] -Indexed[X, 2], Indexed[X, 1];

τ = 0.01;
Block[YY, Y, k1, k2, k3, k4,

 YY = Table[Compile`GetElement[Y, i], i, 1, 2];
 k1 = τ F[YY];
 k2 = τ F[0.5 k1 + YY];
 k3 = τ F[0.5 k2 + YY];
 k4 = τ F[k3 + YY];

 cStep = With[code = YY + (k1 + 2. (k2 + k3) + k4)/6. ,
 Compile[Y, _Real, 1,
 code,
 CompilationTarget -> "C",
 RuntimeOptions -> "Speed"
 ]
 ]
 ];

Now we can apply it 20 million times with NestList and it stills takes only about 2 seconds.

nsteps = 20000000;
xlist = Range[0., τ nsteps, τ];
Ylist = NestList[cStep, 1., 0., nsteps]; // AbsoluteTiming // First

2.08678

Edit

This can be sped up even more my avoiding NestList (the loop behind it can also be compiled which saves several calls to libraries) and by utilizing that the dimension of the ODE is known at compile time. For low dimensional systems, it may be also beneficial to avoid tensor operations altogether and to perform computations in scalar registers as done below.

τ = 0.01;
cFlow = Block[YY, Y, k1, k2, k3, k4, τ, Ylist, j,
 YY = Table[Compile`GetElement[Ylist, j, i], i, 1, 2];
 k1 = τ F[YY];
 k2 = τ F[0.5 k1 + YY];
 k3 = τ F[0.5 k2 + YY];
 k4 = τ F[k3 + YY];
 With[
 code1 = (YY + (k1 + 2. (k2 + k3) + k4)/6)[[1]],
 code2 = (YY + (k1 + 2. (k2 + k3) + k4)/6)[[2]]
 ,
 Compile[Y0, _Real, 1, τ, _Real, n, _Integer,
 Block[Ylist,
 Ylist = Table[0., n + 1, Length[Y0]];
 Ylist[[1]] = Y0;
 Do[
 Ylist[[j + 1, 1]] = code1;
 Ylist[[j + 1, 2]] = code2;
 ,
 j, 1, n];
 Ylist
 ],
 CompilationTarget -> "C", RuntimeOptions -> "Speed"
 ]
 ]
 ];
Ylist2 = cFlow[1., 0., τ, nsteps]; // AbsoluteTiming // First

1.06549

Don't be too upset by parts of the code being highlighted in red; this is on purpose.