Solving a system of ODEs with the Runge-Kutta method

Here is a functional approach. The following will give you one step of the Runge-Kutta formula:

RungeKutta[func_List, yinit_List, y_List, step_] := 
 Module[{k1, k2, k3, k4},
  k1 = step N[func /. MapThread[Rule, {y, yinit}]];
  k2 = step N[func /. MapThread[Rule, {y, k1/2 + yinit}]];
  k3 = step N[func /. MapThread[Rule, {y, k2/2 + yinit}]];
  k4 = step N[func /. MapThread[Rule, {y, k3 + yinit}]];
  yinit + Total[{k1, 2 k2, 2 k3, k4}]/6]

Here, func is a list of functions, yinit a list of initial values, y a list of function variables (in your case that will be {x, y, m, k, s}, and step is the step size of the numerical simulation. You can then use something like NestList to iterate it many times like this:

NestList[RungeKutta[func, #, y, step] &, N[yinit], Round[t/step]]

Here, t is the maximum value of your independent variable. You can also include conditionals to check that the lists are of equal length.

andre has linked you to the method plug-in framework in the comments, but there is a more direct way to implement classical Runge-Kutta, just by supplying its Butcher tableau to "ExplicitRungeKutta". Here's an adaptation from the docs:

ClassicalRungeKuttaCoefficients[4, prec_] := With[{amat = {{1/2}, {0, 1/2}, {0, 0, 1}},
           bvec = {1/6, 1/3, 1/3, 1/6}, cvec = {1/2, 1/2, 1}}, N[{amat, bvec, cvec}, prec]]

{xf, yf} = {x, y} /. First @ 
           NDSolve[{x'[t] == -y[t], y'[t] == x[t], x[0] == 1, y[0] == 0},
                   {x, y}, {t, 0, 6}, 
                   Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4,
                              "Coefficients" -> ClassicalRungeKuttaCoefficients},
                   StartingStepSize -> 1/2];

To obtain the values computed by RK, you can then do this:

xl = MapThread[Append, {xf["Grid"], xf["ValuesOnGrid"]}]
   {{0., 1.}, {0.5, 0.877604}, {1., 0.540588}, {1.5, 0.0714256}, {2., -0.415108},
    {2.5, -0.800012}, {3., -0.989166}, {3.5, -0.936349}, {4., -0.65453}, {4.5, -0.212684},
    {5., 0.281088}, {5.5, 0.706007}, {6., 0.95816}}

yl = MapThread[Append, {yf["Grid"], yf["ValuesOnGrid"]}]
   {{0., 0.}, {0.5, 0.479167}, {1., 0.841037}, {1.5, 0.99713}, {2., 0.90931},
    {2.5, 0.599108}, {3., 0.142441}, {3.5, -0.348969}, {4., -0.754924}, {4.5, -0.976153},
    {5., -0.958587}, {5.5, -0.706572}, {6., -0.281796}}

Plot them:

ListLinePlot[{xl, yl}, Mesh -> All, MeshStyle -> PointSize[Medium]]