Find y(0.4) for `y''=xz^2-y^2`, `x_0=0, y_0=1, z_0=0`, with step length 0.2 using Runge-Kutta 4 method (second order differential equation)

Solution:
Given `y^('')=xz^2-y^2, y(0)=1, y'(0)=0, h=0.2, y(0.4)=?`

put `(dy)/(dx)=z` and differentiate w.r.t. x, we obtain `(d^2y)/(dx^2)=(dz)/(dx)`

We have system of equations
`(dy)/(dx)=z=f(x,y,z)`

`(dz)/(dx)=xz^2-y^2=g(x,y,z)`

Fourth order Runge-Kutta (RK4) method for second order differential equation formula
`k_1=f(x_n,y_n,z_n)`

`l_1=g(x_n,y_n,z_n)`

`k_2=f(x_n+h/2,y_n+(hk_1)/2,z_n+(hl_1)/2)`

`l_2=g(x_n+h/2,y_n+(hk_1)/2,z_n+(hl_1)/2)`

`k_3=f(x_n+h/2,y_n+(hk_2)/2,z_n+(hl_2)/2)`

`l_3=g(x_n+h/2,y_n+(hk_2)/2,z_n+(hl_2)/2)`

`k_4=f(x_n+h,y_n+hk_3,z_n+hl_3)`

`l_4=g(x_n+h,y_n+hk_3,z_n+hl_3)`

`y_(n+1)=y_n+h/6(k_1+2k_2+2k_3+k_4)`

`z_(n+1)=z_n+h/6(l_1+2l_2+2l_3+l_4)`

---

for `n=0,x_0=0,y_0=1,z_0=0`

`k_1=f(x_0,y_0,z_0)`

`=f(0,1,0)`

`=0`

`l_1=g(x_0,y_0,z_0)`

`=g(0,1,0)`

`=-1`

`k_2=f(x_0+h/2,y_0+(hk_1)/2,z_0+(hl_1)/2)`

`=f(0.1,1,-0.1)`

`=-0.1`

`l_2=g(x_0+h/2,y_0+(hk_1)/2,z_0+(hl_1)/2)`

`=g(0.1,1,-0.1)`

`=-0.999`

`k_3=f(x_0+h/2,y_0+(hk_2)/2,z_0+(hl_2)/2)`

`=f(0.1,0.99,-0.0999)`

`=-0.0999`

`l_3=g(x_0+h/2,y_0+(hk_2)/2,z_0+(hl_2)/2)`

`=g(0.1,0.99,-0.0999)`

`=-0.9791`

`k_4=f(x_0+h,y_0+hk_3,z_0+hl_3)`

`=f(0.2,0.98,-0.1958)`

`=-0.1958`

`l_4=g(x_0+h,y_0+hk_3,z_0+hl_3)`

`=g(0.2,0.98,-0.1958)`

`=-0.9528`

Now,
`y_1=y_0+h/6(k_1+2k_2+2k_3+k_4)`

`=1+0.2/6[0+2(-0.1)+2(-0.0999)+(-0.1958)]`

`=0.9801`

`z_1=z_0+h/6(l_1+2l_2+2l_3+l_4)`

`=0+0.2/6[-1+2(-0.999)+2(-0.9791)+(-0.9528)]`

`=-0.197`

`x_1=x_0+h=0+0.2=0.2`

---

for `n=1,x_1=0.2,y_1=0.9801,z_1=-0.197`

`k_1=f(x_1,y_1,z_1)`

`=f(0.2,0.9801,-0.197)`

`=-0.197`

`l_1=g(x_1,y_1,z_1)`

`=g(0.2,0.9801,-0.197)`

`=-0.9529`

`k_2=f(x_1+h/2,y_1+(hk_1)/2,z_1+(hl_1)/2)`

`=f(0.3,0.9604,-0.2923)`

`=-0.2923`

`l_2=g(x_1+h/2,y_1+(hk_1)/2,z_1+(hl_1)/2)`

`=g(0.3,0.9604,-0.2923)`

`=-0.8968`

`k_3=f(x_1+h/2,y_1+(hk_2)/2,z_1+(hl_2)/2)`

`=f(0.3,0.9509,-0.2866)`

`=-0.2866`

`l_3=g(x_1+h/2,y_1+(hk_2)/2,z_1+(hl_2)/2)`

`=g(0.3,0.9509,-0.2866)`

`=-0.8796`

`k_4=f(x_1+h,y_1+hk_3,z_1+hl_3)`

`=f(0.4,0.9228,-0.3729)`

`=-0.3729`

`l_4=g(x_1+h,y_1+hk_3,z_1+hl_3)`

`=g(0.4,0.9228,-0.3729)`

`=-0.796`

Now,
`y_2=y_1+h/6(k_1+2k_2+2k_3+k_4)`

`=0.9801+0.2/6[-0.197+2(-0.2923)+2(-0.2866)+(-0.3729)]`

`=0.9226`

`x_2=x_1+h=0.2+0.2=0.4`

`:.y(0.4)=0.9226`

`n`	`x_n`	`y_n`	`z_n`	`k_1`	`l_1`	`k_2`	`l_2`	`k_3`	`l_3`	`k_4`	`l_4`	`x_(n+1)`	`y_(n+1)`	`z_(n+1)`
0	0	1	0	0	-1	-0.1	-0.999	-0.0999	-0.9791	-0.1958	-0.9528	0.2	0.9801	-0.197
1	0.2	0.9801	-0.197	-0.197	-0.9529	-0.2923	-0.8968	-0.2866	-0.8796	-0.3729	-0.796	0.4	0.9226

---

---