# Solving ODE with Runge-Kutta method [duplicate]

I am trying to solve y''+4y'+3y=0; y[0]=1, y'[0]=-1 with Runge-Kutta method and I have a problem that its DSolve gives the result:

s = DSolve[{y''[x] + 4 y'[x] + 3 y[x] == 0, y[0] == 1, y'[0] == -1},  y[x], x]

Plot[Evaluate[y[x] /. s], {x, 0, 10}, PlotStyle -> Automatic]


And my code gives a very strange not correct result.

n = 100;
h = 10/n;

y = Table[0, n + 1]; z = Table[0, n + 1];
xxx = Table[0, n + 1];

xxx[[1]] = 0;   y[[1]] = 1;  z[[1]] = -1;

fy[xe_, ye_, ze_] = ze;
fz[x_, yt_, zt_] = -3 yt - 4 zt;

For[i = 1, i <= n, i++,

Clear[k1z, k2z, k3z, k1y, k2y, k3y];
yy = Solve[k1y == fy[x + h*1/3, y[[i]] + 1/3*h*k1y, z[[i]]] &&
k2y == fy[x + h*1, y[[i]] + 1/3*h*k1y + 2/3*h*k2y, z[[i]]] &&
k3y == fy[x + h*1, y[[i]] + 0 + 0 + 1*h*k3y, z[[i]]], {k1y, k2y,
k3y}];

k1y0 = k1y /. yy[[1]]; k2y0 = k2y /. yy[[1]]; k3y0 = k3y /. yy[[1]];

zz = Solve[k1z == fy[x + h*1/3, y[[i]], z[[i]] + 1/3*h*k1z] &&
k2z == fy[x + h*1, y[[i]], z[[i]] + 1/3*h*k1z + 2/3*h*k2z] &&
k3z == fy[x + h*1, y[[i]], z[[i]] + 0 + 0 + 1*h*k3z], {k1z, k2z,
k3z}];

k1z0 = k1z /. zz[[1]]; k2z0 = k2z /. zz[[1]]; k3z0 = k3z /. zz[[1]];

y[[i + 1]] = y[[i]] + h (3/4 k1y0 + 3/4 k2y0 - 1/2 k3y0);
z[[i + 1]] = z[[i]] + h (3/4 k1z0 + 3/4 k2z0 - 1/2 k3z0);

xxx[[i + 1]] = xxx[[i]] + h; ]

Print[ListPlot[Transpose[{xxx, y}]]]


I have in the result:

So obviously my code is wrong, can anybody help me why?

Maybe somebody can apply new Runge-Kutta or another numerical method to solve it, I would be very happy to see, how it is working?

• "[A]nother numerical method to solve it": Try NDSolve[]. – Michael E2 May 13 '20 at 18:01
• @Szabolcs unfortunately no, because there is first order of DE, and I need RK for the second order. I think that the problem is somewhere here when I apply RK method for two variables (y, y') – Yura Holubeu May 13 '20 at 18:59
• @MichaelE2 no, this is not an option for me – Yura Holubeu May 13 '20 at 18:59
• My answer there shows a solution for a 2nd order equation (harmonic oscillator). – Szabolcs May 13 '20 at 19:35
• Then what do you mean by "another numerical method to solve it"? – Michael E2 May 13 '20 at 19:41

This code is original, but it can't be used for nonlinear case. Also it is not looks like RK3. Nevertheless after correcting few typos we have

n = 150;
h = 5/n;

y = Table[0, n + 1]; z = Table[0, n + 1];
xxx = Table[0, n + 1];

xxx[[1]] = 0; y[[1]] = 1; z[[1]] = -1;

fy[xe_, ye_, ze_] := ze;
fz[x_, yt_, zt_] := -3 yt - 4 zt;

Do[
yy = Solve[
k1y == fy[xxx[[i]] + h*1/3, y[[i]] + 1/3*h*k1y, z[[i]]] &&
k2y == fy[xxx[[i]] + h*1, y[[i]] + 1/3*h*k1y + 2/3*h*k2y,
z[[i]]] &&
k3y == fy[xxx[[i]] + h*1, y[[i]] + 0 + 0 + 1*h*k3y, z[[i]]], {k1y,
k2y, k3y}];
k1y0 = k1y /. yy[[1]]; k2y0 = k2y /. yy[[1]]; k3y0 = k3y /. yy[[1]];
zz = Solve[
k1z == fz[xxx[[i]] + h*1/3, y[[i]], z[[i]] + 1/3*h*k1z] &&
k2z == fz[xxx[[i]] + h*1, y[[i]],
z[[i]] + 1/3*h*k1z + 2/3*h*k2z] &&
k3z == fz[xxx[[i]] + h*1, y[[i]], z[[i]] + 0 + 0 + 1*h*k3z], {k1z,
k2z, k3z}];
k1z0 = k1z /. zz[[1]]; k2z0 = k2z /. zz[[1]]; k3z0 = k3z /. zz[[1]];
y[[i + 1]] = y[[i]] + h (3/4 k1y0 + 3/4 k2y0 - 1/2 k3y0);
z[[i + 1]] = z[[i]] + h (3/4 k1z0 + 3/4 k2z0 - 1/2 k3z0);
xxx[[i + 1]] = xxx[[i]] + h;, {i, 1, n}]


Comparison with exact solution $$e^{-x}$$

Show[Plot[Exp[-x], {x, 0, 5}, PlotRange -> All],
ListPlot[Transpose[{xxx, y}], PlotRange -> All, PlotStyle -> Red]]


• woow, it seems you are genius! Can you please tell me where was my mistake in my code? – Yura Holubeu May 13 '20 at 20:58
• can you please tell me, if instead of fz[x_, yt_, zt_] := -3 yt - 4 zt; I put in your code fz[x_, yt_, zt_] := 11yt +10 zt; everything falls and solution begans to grow? I just have no idea.... – Yura Holubeu May 13 '20 at 21:06
• @YuraHolubeu Actually this is your code, not mine. I just made it working. You use x instead of xxx[[i]] and fy instead of fz. Do you try to implement RK3 schema? – Alex Trounev May 13 '20 at 22:20
• yes, it was so stupid that I used x, I couldn't find this for hours... But also I am trying to figure out why I had such big numbers?... because frankly speaking my original problem was with fz[x_, yt_, zt_] := 11yt +10 zt; and RK should give Exp[-x]... I have no idea, why I have such a sharp rise... – Yura Holubeu May 13 '20 at 22:46
• and unfortunately if I put fz[x_, yt_, zt_] := x Sqrt[yt];  , then everything begins to count for a really long time, while it should take seconds. I have no idea what is wrong with RK???? Maybe you could know? because probabity that I will figure out it on my own is around zero... – Yura Holubeu May 13 '20 at 23:45