Question: How do I solve ODE system in RK method?

How to solve and plot a ODE system in RK method.
eq1 := diff(f(x), x, x, x)-(1/2)*Sc*sin(alpha)*g(x)*(diff(g(x), x, x))+(1/2)*x*cos(alpha)*(diff(f(x), x, x))+(1/2)*sin(alpha)*f(x)*(diff(f(x), x, x)) = 0; eq2 := (diff(g(x), x, x, x))/Pm+(1/2)*x*cos(alpha)*(diff(g(x), x, x))+sin(alpha)*f(x)*(diff(g(x), x, x))-sin(alpha)*(diff(f(x), x, x))*g(x) = 0; eq3 := (diff(theta(x), x, x))/Pr+(1/2)*x*cos(alpha)*(diff(theta(x), x))+(1/2)*x*(diff(f(x), x))*(diff(theta(x), x))+sin(alpha)*(x*(diff(f(x), x))-f(x))*(diff(theta(x), x))-Nb*(diff(s(x), x))*(diff(theta(x), x))-Nt*(diff(theta(x), x))^2+(1/4)*Sc*Br*sin(alpha)^2*(diff(f(x), x))^2*(x*(diff(g(x), x))-g(x))+(diff(g(x), x))^2*(x*(diff(f(x), x))-f(x)) = 0; eq4 := diff(s(x), x, x)+S*((1/2)*cos(alpha)*x*(diff(s(x), x))+(1/2)*sin(alpha)*f(x)*(diff(s(x), x)))+Nt*(diff(theta(x), x, x))/Nb = 0

ics := f(0) = 0, (D(f))(0) = 1, g(0) = 0, (D(g))(0) = 1, theta(0) = 1, s(0) = 1; bcs := (D(f))(100) = 0, (D(g))(100) = 0, theta(100) = 0, s(100) = 0

alpha = - 30 degree, Sc = 1.0, Pm = .1, Pr = 6.2, Nb = .1, Nt = .1, Br = .5, S = 1