Re your point [14]: Using f as argument to RK4step and RK4 would have two advantages:

1. You would not have to use the name f, which is currently hard-coded into RK4step and RK4, for your function prescribing the ODE's right-hand side. You could call it anything you like!
2. It would make these function more like DSolve and NDSolve that are built into *Mathematica": in those two, the first argument is always the differential equation (or list consisting of differential equation and initial or boundary condtions) to be solved.

This would still not be quite like the built-in functions, in that the first argument would be just the name of a function (or, alternatively, a "pure function" expressed in a Function[....] expression or the equivalent abbreviation using # and &).

To be much more like the built-in functions, you would need to allow the first argument to be an entire "function expression" including variables, e.g., (in the case of a scalar ODE, for simplicity):

x^2 - y

But this makes coding RK4step and RK4 trickier, for you would also then need a new 2nd argument that is a list, in this example, {x, y}, of the names of the variables in that first argument.

In general, it looks like you are treating a Mathematica notebook, and the Wolfram Language code used in it, as if it were a Matlab .m file and Matlab code. In contrast to a Matlab .m file, one does not ordinarily evaluate the entire notebook, but rather evaluate cells one-by-one, or in groups, as needed.

In your notebook, you are calling RK4 within the For loop before you define the function RK4. (Note also that RK4 calls the function f before you define f.) That is why you need to evaluate the entire notebook twice — unless, that is, you evaluate cells individually but in the order of defining f, defining RK4, and only then the For loop.

Aside from that, there are some issues with your code.

• There is utterly no need for a For loop to print what you want. You should be able to do this with Table or similar constructs.
• The function RK4 does not return a result as output (well, technically, it returns Null because the last expression within it ends with a semicolon.
• Similarly, in the definition of f, you have no actual output but rather perniciously assign values to the (first two positions of) the global y.
• In the definition of f, you have a local variable i which is not subsequently used, so why have it?
• At several places in the definition of RK4 you simply invoke f at various input pairs but do nothing with them. For example, immediately after your list of local variables in Module, you have f[x, yn];. This evaluates f at the current values of x and yn and then in effect just throws that value of f away, because you don't assign it to anything; perhaps you really meant y = f[x, yn] there. And then you could simply set k1 = h f[x, yn].
• This is a matter of taste, perhaps, but there is utterly no need to insert the symbol * to indicate multiplication; instead, just type a space between, e.g., h and y (so as not to confuse the product with a new symbol hy); Mathematica will provide an appropriate thin space there so it looks nice!
• Your function RK4 takes as inputs an x , a yn, an m, and an h. Why, then, do you bother to define them globally at the beginning?

• Inside RK4, the way you create the coefficients k1, k2, k3, k4 is unnecessarily verbose. This can be done much more succinctly — and clearly! — as follows:

  k1 = k2 = k3 = k4 = ConstantArray[0., m]