Solving ODE with Runge-Kutta method [duplicate]

Question

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?

Answer 1

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]]