user
Given the example Differential equation:

:   y
    ′
    (
    t
    )
    =
    t
    ⋅


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

With initial condition:

:   t

    0
    =
    0
    {\displaystyle {\displaystyle t\_{0}=0}}
     and 

    y

    0
    =
    y
    (

    t

    0
    )
    =
    y
    (
    0
    )
    =
    1
    {\displaystyle {\displaystyle y\_{0}=y(t\_{0})=y(0)=1}}

This equation has an exact solution:

:   y
    (
    t
    )
    =

    1
    16
    (

    t

    2
    +
    4

    )

    2
    {\displaystyle {\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 {\displaystyle t=0\ldots 10}}
   with a step value of 

  δ
  t
  =
  0.1
  {\displaystyle {\displaystyle \delta t=0.1}}
   (101 total points, the first being given)
* Print the calculated values of 

  y
  {\displaystyle {\displaystyle y}}
   at whole numbered 

  t
  {\displaystyle {\displaystyle t}}
  's (

  0.0
  ,
  1.0
  ,
  …
  10.0
  {\displaystyle {\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 {\displaystyle y\_{n}}}
 and 

t

n
{\displaystyle {\displaystyle t\_{n}}}
 calculate:

:   δ

    y

    1
    =
    δ
    t
    ⋅

    y
    ′
    (

    t

    n
    ,

    y

    n
    )

    {\displaystyle {\displaystyle \delta y\_{1}=\delta t\cdot y'(t\_{n},y\_{n})\quad }}
:   δ

    y

    2
    =
    δ
    t
    ⋅

    y
    ′
    (

    t

    n
    +

    1
    2
    δ
    t
    ,

    y

    n
    +

    1
    2
    δ

    y

    1
    )
    {\displaystyle {\displaystyle \delta y\_{2}=\delta t\cdot 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 {\displaystyle \delta y\_{3}=\delta t\cdot 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 {\displaystyle \delta y\_{4}=\delta t\cdot 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 {\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 {\displaystyle t\_{n+1}=t\_{n}+\delta t\quad }}
assistant
(formerly Perl 6)

**Works with**: rakudo version 2016.03

```
sub runge-kutta(&yp) {
    return -> \t, \y, \δt {
        my $a = δt * yp( t, y );
        my $b = δt * yp( t + δt/2, y + $a/2 );
        my $c = δt * yp( t + δt/2, y + $b/2 );
        my $d = δt * yp( t + δt, y + $c );
        ($a + 2*($b + $c) + $d) / 6;
    }
}
 
constant δt = .1;
my &δy = runge-kutta { $^t * sqrt($^y) };
 
loop (
    my ($t, $y) = (0, 1);
    $t <= 10;
    ($t, $y) »+=« (δt, δy($t, $y, δt))
) {
    printf "y(%2d) = %12f ± %e\n", $t, $y, abs($y - ($t**2 + 4)**2 / 16)
    if $t %% 1;
}

```

Output:

```
y( 0) =     1.000000 ± 0.000000e+00
y( 1) =     1.562500 ± 1.457219e-07
y( 2) =     3.999999 ± 9.194792e-07
y( 3) =    10.562497 ± 2.909562e-06
y( 4) =    24.999994 ± 6.234909e-06
y( 5) =    52.562489 ± 1.081970e-05
y( 6) =    99.999983 ± 1.659460e-05
y( 7) =   175.562476 ± 2.351773e-05
y( 8) =   288.999968 ± 3.156520e-05
y( 9) =   451.562459 ± 4.072316e-05
y(10) =   675.999949 ± 5.098329e-05
```