Given the example Differential equation:






y
′

(
t
)
=
t
×


y
(
t
)




{\displaystyle y'(t)=t\times {\sqrt {y(t)}}}


With initial condition:






t

0


=
0


{\displaystyle t_{0}=0}

 and 




y

0


=
y
(

t

0


)
=
y
(
0
)
=
1


{\displaystyle y_{0}=y(t_{0})=y(0)=1}


This equation has an exact solution:





y
(
t
)
=



1
16



(

t

2


+
4

)

2




{\displaystyle y(t)={\tfrac {1}{16}}(t^{2}+4)^{2}}




Task
Demonstrate the commonly used explicit   fourth-order Runge–Kutta method   to solve the above differential equation.

 Solve the given differential equation over the range 



t
=
0
…
10


{\displaystyle t=0\ldots 10}

 with a step value of 



δ
t
=
0.1


{\displaystyle \delta t=0.1}

 (101 total points, the first being given)
 Print the calculated values of 



y


{\displaystyle y}

 at whole numbered 



t


{\displaystyle t}

's (



0.0
,
1.0
,
…
10.0


{\displaystyle 0.0,1.0,\ldots 10.0}

) along with error as compared to the exact solution.


Method summary
Starting with a given 




y

n




{\displaystyle y_{n}}

 and 




t

n




{\displaystyle t_{n}}

 calculate:





δ

y

1


=
δ
t
×

y
′

(

t

n


,

y

n


)



{\displaystyle \delta y_{1}=\delta t\times y'(t_{n},y_{n})\quad }






δ

y

2


=
δ
t
×

y
′

(

t

n


+



1
2



δ
t
,

y

n


+



1
2



δ

y

1


)


{\displaystyle \delta y_{2}=\delta t\times y'(t_{n}+{\tfrac {1}{2}}\delta t,y_{n}+{\tfrac {1}{2}}\delta y_{1})}






δ

y

3


=
δ
t
×

y
′

(

t

n


+



1
2



δ
t
,

y

n


+



1
2



δ

y

2


)


{\displaystyle \delta y_{3}=\delta t\times y'(t_{n}+{\tfrac {1}{2}}\delta t,y_{n}+{\tfrac {1}{2}}\delta y_{2})}






δ

y

4


=
δ
t
×

y
′

(

t

n


+
δ
t
,

y

n


+
δ

y

3


)



{\displaystyle \delta y_{4}=\delta t\times y'(t_{n}+\delta t,y_{n}+\delta y_{3})\quad }


then:






y

n
+
1


=

y

n


+



1
6



(
δ

y

1


+
2
δ

y

2


+
2
δ

y

3


+
δ

y

4


)


{\displaystyle y_{n+1}=y_{n}+{\tfrac {1}{6}}(\delta y_{1}+2\delta y_{2}+2\delta y_{3}+\delta y_{4})}







t

n
+
1


=

t

n


+
δ
t



{\displaystyle t_{n+1}=t_{n}+\delta t\quad }




function rk4(f, t0, y0, t1, n) {
	h = (t1-t0)/(n-1)
	a = J(n, 2, 0)
	a[1, 1] = t = t0
	a[1, 2] = y = y0
	for (i=2; i<=n; i++) {
		k1 = h*(*f)(t, y)
		k2 = h*(*f)(t+0.5*h, y+0.5*k1)
		k3 = h*(*f)(t+0.5*h, y+0.5*k2)
		k4 = h*(*f)(t+h, y+k3)
		t = t+h
		y = y+(k1+2*k2+2*k3+k4)/6
		a[i, 1] = t
		a[i, 2] = y
	}
	return(a)
}
 
function f(t, y) {
	return(t*sqrt(y))
}
 
a = rk4(&f(), 0, 1, 10, 101)
t = a[., 1]
a = a, a[., 2]:-(t:^2:+4):^2:/16
a[range(1,101,10), .]
 
                   1              2              3
     +----------------------------------------------+
   1 |             0              1              0  |
   2 |             1    1.562499854   -1.45722e-07  |
   3 |             2    3.999999081   -9.19479e-07  |
   4 |             3    10.56249709   -2.90956e-06  |
   5 |             4    24.99999377   -6.23491e-06  |
   6 |             5    52.56248918   -.0000108197  |
   7 |             6    99.99998341   -.0000165946  |
   8 |             7    175.5624765   -.0000235177  |
   9 |             8    288.9999684   -.0000315652  |
  10 |             9    451.5624593   -.0000407232  |
  11 |            10     675.999949   -.0000509833  |
     +----------------------------------------------+