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 }




begin
    real procedure rk4 ( real procedure f ; real value y, x, dx ) ;
    begin  % Fourth-order Runge-Kutta method %
       real dy1, dy2, dy3, dy4;
       dy1 := dx * f(x, y);
       dy2 := dx * f(x + dx / 2.0, y + dy1 / 2.0);
       dy3 := dx * f(x + dx / 2.0, y + dy2 / 2.0);
       dy4 := dx * f(x + dx, y + dy3);
       y + (dy1 + 2.0 * dy2 + 2.0 * dy3 + dy4) / 6.0
    end rk4;
    real x0, x1, y0, dx;
    integer numPoints;
    x0 := 0; x1 := 10; y0 := 1.0;                                % Boundary conditions. %
    dx := 0.1;                                                   % Step size. %
    numPoints := entier ((x1 - x0) / dx + 0.5);                  % Add 0.5 for rounding errors. %
    begin
        real procedure dyByDx ( real value x, y ) ; x * sqrt(y); % Differential equation. %
        real array y ( 0 :: numPoints); y(0) := y0;              % Grid and starting point. %
        for i := 1 until numPoints do y(i) := rk4 (dyByDx, y(i-1), x0 + dx * (i - 1), dx);
        write( "   x              true y         calc y       relative error" );
        for i := 0 step 10 until numPoints do begin
            real x, trueY;
            x     := x0 + dx * i;
            trueY := (x * x + 4.0) ** 2 / 16.0;
            write( r_format := "A", r_w := 12, r_d := 7, s_w := 3, x, trueY, y( i )
                 , r_format := "S", r_w := 12, y( i ) / trueY - 1
                 )
        end for_i
    end
end.