(* r.h.s. function f, step - size h, intial values for x, y *)
f[x_, y_] := x y^(1/3); h = 0.1; x = 1; y = 1;

(* The classical RK4 scheme *)
Do[
 Print[{x, y}];
 k1 = h f[x, y]; k2 = h f[x + 0.5 h, y + 0.5 k1]; 
 k3 = h f[x + 0.5 h, y + 0.5 k2]; 
 k4 = h f[x + 1 h, y + k3]; (* observe that h is a factor when computing the k-values, this is the classical notation *)
 Print[{k1, k2, k3, k4}];
 x = x + h;
 y = y + (k1 + 2 k2 + 2 k3 + k4)/6; (* must be replaced by the exact y-value when computing local errors ! *)
 ,
 5
]