Solve the following differential equation using the fourth-order Runge-Kutta method [duplicate]

I need to use Mathematica to solve the following DE by the fourth-order Runge-Kutta method and find $y(10)$: $$y'(x)=\frac{1}{2y}, y(0)=1$$ The following is the code I have so far (with step size h=0.5):

Clear[x, y, n, h, S1, S2, S3, S4, i];
y[0] = 1; y[1] = Sqrt[2]; h = 0.5; n = 20;
Do[x[i] = 0 + (i - 1)*h, {i, 1, n + 1}]

y'[x_] = f[x_, y_] = 1/(2*y);
Do[{S1 = f[x[i], y[i]],
S2 = f[x[i] + h/2, y[i] + h*S1/2],
S3 = f[x[i] + h/2, y[i] + h*S2/2],
S4 = f[x[i] + h, y[i] + h*S3],
y[i + 1] = y[i] + (S1 + 2 S2 + 2 S3 + S4)*h/6}, {i, 1, n}]
Do[Print[x[i], " ", y[i]], {i, 1, n + 1}]

This code results in an approximation of $y(10)=3.4641$, which is fine; however, when I change the step size and the corresponding $n$ value I still get $y(10)=3.4641$ regardless of the value of the step size, which doesn't make sense to me. I'm almost certain that the smaller the step size one chooses, the better the approximation should be, so we shouldn't observe the same value in error. I'm not sure how to proceed from here, but I'm sure there is an error in my code somewhere. Any help is greatly appreciated.

• Why y[1] = Sqrt[2]? May 10, 2017 at 3:05
• initially I used the given initial condition only (y(0) = 1) but Mathematica wouldn't give an output in a timely fashion May 10, 2017 at 3:33

You are suppose to give the initial condition y[1]=1 not y[0]=1 because your loop starts from i=1.

Your this condition y[1]=Sqrt[2] doesn't make any sense.

If you start your loop from i=0 and ends at i=n then you can use y[0]=1?

For h=0.01, I get y[10]=3.31632. Comparing this with the results given by NDSolve,

Clear[x, y];
sol = NDSolve[{y'[x] == 1/(2*y[x]), y[0] == 1}, y[x], {x, 0, 10}];
y[x] /. sol /. x -> 10

3.31662

Clear[x, y, S1, S2];
a = 0; b = 10;
y[0] = 1; h = 0.01; n = (b - a)/h;
Do[x[i] = 0.0 + (i)*h, {i, 0, n}]

f[x_, y_] = 1/(2*y);
Do[{S1 = f[x[i], y[i]],
S2 = f[x[i] + h, y[i] + h*S1],
S3 = f[x[i] + h/2, y[i] + h*S2/2],
S4 = f[x[i] + h, y[i] + h*S3],
y[i + 1] = y[i] + (S1 + 2 S2 + 2 S3 + S4)*h/6}, {i, 0, n}]
Do[Print[x[i], "       ", y[i]], {i, 0, n}];
• can you explain why i must give initial condition y[1]==1 and not the others? May 10, 2017 at 3:56