y = [y, u];
k1 = f(u);
k2 = f(u + k1.*(h/2));
k3 = f(u + k2.*(h/2));
k4 = f(u + k3.*h);

u = u + (k1 + k2.*2 + k3.*2 + k4).*(h/6);