# Solve the following differential equation using the fourth-order Runge-Kutta method

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]? – Michael E2 May 10 '17 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 – linearalgebrathrowaway May 10 '17 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? – linearalgebrathrowaway May 10 '17 at 3:56
• @linearalgebrathrowaway check edit please – zhk May 10 '17 at 3:57