# Runge-Kutta method

Runge-Kutta method
You are encouraged to solve this task according to the task description, using any language you may know.

Given the example Differential equation:

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

With initial condition:

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

This equation has an exact solution:

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

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

• Solve the given differential equation over the range ${\displaystyle t=0\ldots 10}$ with a step value of ${\displaystyle \delta t=0.1}$ (101 total points, the first being given)
• Print the calculated values of ${\displaystyle y}$ at whole numbered ${\displaystyle t}$'s (${\displaystyle 0.0,1.0,\ldots 10.0}$) along with error as compared to the exact solution.
Method summary

Starting with a given ${\displaystyle y_{n}}$ and ${\displaystyle t_{n}}$ calculate:

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

then:

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

  type Floaty is digits 15;
type Floaty_Array is array (Natural range <>) of Floaty;
package FIO is new Ada.Text_IO.Float_IO(Floaty); use FIO;
type Derivative is access function(t, y : Floaty) return Floaty;
package Math is new Ada.Numerics.Generic_Elementary_Functions (Floaty);
function calc_err (t, calc : Floaty) return Floaty;

procedure Runge (yp_func : Derivative; t, y : in out Floaty_Array;
dt : Floaty) is
dy1, dy2, dy3, dy4 : Floaty;
begin
for n in t'First .. t'Last-1 loop
dy1 := dt * yp_func(t(n), y(n));
dy2 := dt * yp_func(t(n) + dt / 2.0, y(n) + dy1 / 2.0);
dy3 := dt * yp_func(t(n) + dt / 2.0, y(n) + dy2 / 2.0);
dy4 := dt * yp_func(t(n) + dt, y(n) + dy3);
t(n+1) := t(n) + dt;
y(n+1) := y(n) + (dy1 + 2.0 * (dy2 + dy3) + dy4) / 6.0;
end loop;
end Runge;

procedure Print (t, y : Floaty_Array; modnum : Positive) is begin
for i in t'Range loop
if i mod modnum = 0 then
Put("y(");   Put (t(i), Exp=>0, Fore=>0, Aft=>1);
Put(") = "); Put (y(i), Exp=>0, Fore=>0, Aft=>8);
Put(" Error:"); Put (calc_err(t(i),y(i)), Aft=>5);
New_Line;
end if;
end loop;
end Print;

  function yprime (t, y : Floaty) return Floaty is begin
return t * Math.Sqrt (y);
end yprime;
function calc_err (t, calc : Floaty) return Floaty is
actual : constant Floaty := (t**2 + 4.0)**2 / 16.0;
begin return abs(actual-calc);
end calc_err;

dt : constant Floaty := 0.10;
N : constant Positive := 100;
t_arr, y_arr : Floaty_Array(0 .. N);


begin

  t_arr(0) := 0.0;
y_arr(0) := 1.0;
Runge (yprime'Access, t_arr, y_arr, dt);
Print (t_arr, y_arr, 10);


end RungeKutta;</lang>

Output:
y(0.0) = 1.00000000 Error: 0.00000E+00
y(1.0) = 1.56249985 Error: 1.45722E-07
y(2.0) = 3.99999908 Error: 9.19479E-07
y(3.0) = 10.56249709 Error: 2.90956E-06
y(4.0) = 24.99999377 Error: 6.23491E-06
y(5.0) = 52.56248918 Error: 1.08197E-05
y(6.0) = 99.99998341 Error: 1.65946E-05
y(7.0) = 175.56247648 Error: 2.35177E-05
y(8.0) = 288.99996843 Error: 3.15652E-05
y(9.0) = 451.56245928 Error: 4.07232E-05
y(10.0) = 675.99994902 Error: 5.09833E-05

## AWK

<lang AWK>

1. syntax: GAWK -f RUNGE-KUTTA_METHOD.AWK
2. converted from BBC BASIC

BEGIN {

   print(" t    y         error")
y = 1
for (i=0; i<=100; i++) {
t = i / 10
if (t == int(t)) {
actual = ((t^2+4)^2) / 16
printf("%2d %12.7f %g\n",t,y,actual-y)
}
k1 = t * sqrt(y)
k2 = (t + 0.05) * sqrt(y + 0.05 * k1)
k3 = (t + 0.05) * sqrt(y + 0.05 * k2)
k4 = (t + 0.10) * sqrt(y + 0.10 * k3)
y += 0.1 * (k1 + 2 * (k2 + k3) + k4) / 6
}
exit(0)


} </lang>

output:

 t    y         error
0    1.0000000 0
1    1.5624999 1.45722e-007
2    3.9999991 9.19479e-007
3   10.5624971 2.90956e-006
4   24.9999938 6.23491e-006
5   52.5624892 1.08197e-005
6   99.9999834 1.65946e-005
7  175.5624765 2.35177e-005
8  288.9999684 3.15652e-005
9  451.5624593 4.07232e-005
10  675.9999490 5.09833e-005


## BASIC

### BBC BASIC

<lang bbcbasic> y = 1.0

     FOR i% = 0 TO 100
t = i% / 10

IF t = INT(t) THEN
actual = ((t^2 + 4)^2) / 16
PRINT "y("; t ") = "; y TAB(20) "Error = ";  actual - y
ENDIF

k1 =  t * SQR(y)
k2 = (t + 0.05) * SQR(y + 0.05 * k1)
k3 = (t + 0.05) * SQR(y + 0.05 * k2)
k4 = (t + 0.10) * SQR(y + 0.10 * k3)
y += 0.1 * (k1 + 2 * (k2 + k3) + k4) / 6
NEXT i%</lang>


Output:

y(0) = 1            Error = 0
y(1) = 1.56249985   Error = 1.45721892E-7
y(2) = 3.99999908   Error = 9.19479201E-7
y(3) = 10.5624971   Error = 2.90956245E-6
y(4) = 24.9999938   Error = 6.23490936E-6
y(5) = 52.5624892   Error = 1.08196974E-5
y(6) = 99.9999834   Error = 1.65945964E-5
y(7) = 175.562476   Error = 2.35177287E-5
y(8) = 288.999968   Error = 3.15652015E-5
y(9) = 451.562459   Error = 4.07231605E-5
y(10) = 675.999949  Error = 5.09832905E-5


## C

<lang c>#include <stdio.h>

1. include <stdlib.h>
2. include <math.h>

double rk4(double(*f)(double, double), double dx, double x, double y) { double k1 = dx * f(x, y), k2 = dx * f(x + dx / 2, y + k1 / 2), k3 = dx * f(x + dx / 2, y + k2 / 2), k4 = dx * f(x + dx, y + k3); return y + (k1 + 2 * k2 + 2 * k3 + k4) / 6; }

double rate(double x, double y) { return x * sqrt(y); }

int main(void) { double *y, x, y2; double x0 = 0, x1 = 10, dx = .1; int i, n = 1 + (x1 - x0)/dx; y = malloc(sizeof(double) * n);

for (y[0] = 1, i = 1; i < n; i++) y[i] = rk4(rate, dx, x0 + dx * (i - 1), y[i-1]);

printf("x\ty\trel. err.\n------------\n"); for (i = 0; i < n; i += 10) { x = x0 + dx * i; y2 = pow(x * x / 4 + 1, 2); printf("%g\t%g\t%g\n", x, y[i], y[i]/y2 - 1); }

return 0;

}</lang>output (errors are relative)
x       y       rel. err.
------------
0       1       0
1       1.5625  -9.3262e-08
2       4       -2.2987e-07
3       10.5625 -2.75462e-07
4       25      -2.49396e-07
5       52.5625 -2.05844e-07
6       100     -1.65946e-07
7       175.562 -1.33956e-07
8       289     -1.09222e-07
9       451.562 -9.01828e-08
10      676     -7.54191e-08


## D

<lang d>import std.stdio, std.math, std.typecons;

alias FP = real; alias FPs = Typedef!(FP[101]);

void runge(in FP function(in FP, in FP)

          pure nothrow @safe @nogc yp_func,
ref FPs t, ref FPs y, in FP dt) pure nothrow @safe @nogc {
foreach (immutable n; 0 .. t.length - 1) {
immutable FP
dy1 = dt * yp_func(t[n], y[n]),
dy2 = dt * yp_func(t[n] + dt / 2.0, y[n] + dy1 / 2.0),
dy3 = dt * yp_func(t[n] + dt / 2.0, y[n] + dy2 / 2.0),
dy4 = dt * yp_func(t[n] + dt, y[n] + dy3);
t[n + 1] = t[n] + dt;
y[n + 1] = y[n] + (dy1 + 2.0 * (dy2 + dy3) + dy4) / 6.0;
}


}

FP calc_err(in FP t, in FP calc) pure nothrow @safe @nogc {

   immutable FP actual = (t ^^ 2 + 4.0) ^^ 2 / 16.0;
return abs(actual - calc);


}

void main() {

   enum FP dt = 0.10;
FPs t_arr, y_arr;

   t_arr[0] = 0.0;
y_arr[0] = 1.0;
runge((t, y) => t * y.sqrt, t_arr, y_arr, dt);

   foreach (immutable i; 0 .. t_arr.length)
if (i % 10 == 0)
writefln("y(%.1f) = %.8f Error: %.6g",
t_arr[i], y_arr[i],
calc_err(t_arr[i], y_arr[i]));


}</lang>

Output:
y(0.0) = 1.00000000 Error: 0
y(1.0) = 1.56249985 Error: 1.45722e-07
y(2.0) = 3.99999908 Error: 9.19479e-07
y(3.0) = 10.56249709 Error: 2.90956e-06
y(4.0) = 24.99999377 Error: 6.23491e-06
y(5.0) = 52.56248918 Error: 1.08197e-05
y(6.0) = 99.99998341 Error: 1.65946e-05
y(7.0) = 175.56247648 Error: 2.35177e-05
y(8.0) = 288.99996843 Error: 3.15652e-05
y(9.0) = 451.56245928 Error: 4.07232e-05
y(10.0) = 675.99994902 Error: 5.09833e-05

## Dart

<lang dart>import 'dart:math' as Math;

num RungeKutta4(Function f, num t, num y, num dt){

 num k1 = dt * f(t,y);
num k2 = dt * f(t+0.5*dt, y + 0.5*k1);
num k3 = dt * f(t+0.5*dt, y + 0.5*k2);
num k4 = dt * f(t + dt, y + k3);
return y + (1/6) * (k1 + 2*k2 + 2*k3 + k4);


}

void main(){

 num t  = 0;
num dt = 0.1;
num tf = 10;
num totalPoints = ((tf-t)/dt).floor()+1;
num y  = 1;
Function f  = (num t, num y) => t * Math.sqrt(y);
Function actual = (num t) => (1/16) * (t*t+4)*(t*t+4);
for (num i = 0; i <= totalPoints; i++){
num relativeError = (actual(t) - y)/actual(t);
if (i%10 == 0){
print('y(${t.round().toStringAsPrecision(3)}) =${y.toStringAsPrecision(11)}  Error = ${relativeError.toStringAsPrecision(11)}'); } y = RungeKutta4(f, t, y, dt); t += dt; }  }</lang> Output: y(0.00) = 1.0000000000 Error = 0.0000000000 y(1.00) = 1.5624998543 Error = 9.3262010950e-8 y(2.00) = 3.9999990805 Error = 2.2986980086e-7 y(3.00) = 10.562497090 Error = 2.7546153479e-7 y(4.00) = 24.999993765 Error = 2.4939637555e-7 y(5.00) = 52.562489180 Error = 2.0584442034e-7 y(6.00) = 99.999983405 Error = 1.6594596090e-7 y(7.00) = 175.56247648 Error = 1.3395644308e-7 y(8.00) = 288.99996843 Error = 1.0922214534e-7 y(9.00) = 451.56245928 Error = 9.0182772312e-8 y(10.0) = 675.99994902 Error = 7.5419063100e-8  ## Fortran <lang fortran>program rungekutta implicit none real(kind=kind(1.0D0)) :: t,dt,tstart,tstop real(kind=kind(1.0D0)) :: y,k1,k2,k3,k4 tstart =0.0D0 ; tstop =10.0D0 ; dt = 0.1D0 y = 1.0D0 t = tstart write(6,'(A,f4.1,A,f12.8,A,es13.6)') 'y(',t,') = ',y,' Error = '&  &,abs(y-(t**2+4.0d0)**2/16.0d0)  do; if ( t .ge. tstop ) exit  k1 = f (t , y ) k2 = f (t+0.5D0 * dt, y +0.5D0 * dt * k1) k3 = f (t+0.5D0 * dt, y +0.5D0 * dt * k2) k4 = f (t+ dt, y + dt * k3) y = y + dt *( k1 + 2.0D0 *( k2 + k3 ) + k4 )/6.0D0 t = t + dt if(abs(real(nint(t))-t) .le. 1.0D-12) then write(6,'(A,f4.1,A,f12.8,A,es13.6)') 'y(',t,') = ',y,' Error = '& &,abs(y-(t**2+4.0d0)**2/16.0d0) end if  end do contains  function f (t,y) implicit none real(kind=kind(1.0D0)),intent(in) :: y,t real(kind=kind(1.0D0)) :: f f = t*sqrt(y) end function f  end program rungekutta </lang> Output: y( 0.0) = 1.00000000 Error = 0.000000E+00 y( 1.0) = 1.56249985 Error = 1.457219E-07 y( 2.0) = 3.99999908 Error = 9.194792E-07 y( 3.0) = 10.56249709 Error = 2.909562E-06 y( 4.0) = 24.99999377 Error = 6.234909E-06 y( 5.0) = 52.56248918 Error = 1.081970E-05 y( 6.0) = 99.99998341 Error = 1.659460E-05 y( 7.0) = 175.56247648 Error = 2.351773E-05 y( 8.0) = 288.99996843 Error = 3.156520E-05 y( 9.0) = 451.56245928 Error = 4.072316E-05 y(10.0) = 675.99994902 Error = 5.098329E-05  ## F# <lang F Sharp> open System let y'(t,y) = t * sqrt(y) let RungeKutta4 t0 y0 t_max dt =  let dy1(t,y) = dt * y'(t,y) let dy2(t,y) = dt * y'(t+dt/2.0, y+dy1(t,y)/2.0) let dy3(t,y) = dt * y'(t+dt/2.0, y+dy2(t,y)/2.0) let dy4(t,y) = dt * y'(t+dt, y+dy3(t,y))   (t0,y0) |> Seq.unfold (fun (t,y) -> if ( t <= t_max) then Some((t,y), (Math.Round(t+dt, 6), y + ( dy1(t,y) + 2.0*dy2(t,y) + 2.0*dy3(t,y) + dy4(t,y))/6.0)) else None )  let y_exact t = (pown (pown t 2 + 4.0) 2)/16.0 RungeKutta4 0.0 1.0 10.0 0.1  |> Seq.filter (fun (t,y) -> t % 1.0 = 0.0 ) |> Seq.iter (fun (t,y) -> Console.WriteLine("y({0})={1}\t(relative error:{2})", t, y, (y / y_exact(t))-1.0) )</lang>  Output: y(0)=1 (relative error:0) y(1)=1.56249985427811 (relative error:-9.32620110027926E-08) y(2)=3.9999990805208 (relative error:-2.29869800194571E-07) y(3)=10.5624970904376 (relative error:-2.75461533583155E-07) y(4)=24.9999937650906 (relative error:-2.49396374552013E-07) y(5)=52.5624891803026 (relative error:-2.05844421730106E-07) y(6)=99.9999834054036 (relative error:-1.65945964192282E-07) y(7)=175.562476482271 (relative error:-1.33956447156969E-07) y(8)=288.999968434799 (relative error:-1.09222150213029E-07) y(9)=451.56245927684 (relative error:-9.01827772459285E-08) y(10)=675.99994901671 (relative error:-7.54190684348899E-08)  ## Go Works with: Go1 <lang go>package main import (  "fmt" "math"  ) type ypFunc func(t, y float64) float64 type ypStepFunc func(t, y, dt float64) float64 // newRKStep takes a function representing a differential equation // and returns a function that performs a single step of the forth-order // Runge-Kutta method. func newRK4Step(yp ypFunc) ypStepFunc {  return func(t, y, dt float64) float64 { dy1 := dt * yp(t, y) dy2 := dt * yp(t+dt/2, y+dy1/2) dy3 := dt * yp(t+dt/2, y+dy2/2) dy4 := dt * yp(t+dt, y+dy3) return y + (dy1+2*(dy2+dy3)+dy4)/6 }  } // example differential equation func yprime(t, y float64) float64 {  return t * math.Sqrt(y)  } // exact solution of example func actual(t float64) float64 {  t = t*t + 4 return t * t / 16  } func main() {  t0, tFinal := 0, 10 // task specifies times as integers, dtPrint := 1 // and to print at whole numbers. y0 := 1. // initial y. dtStep := .1 // step value.   t, y := float64(t0), y0 ypStep := newRK4Step(yprime) for t1 := t0 + dtPrint; t1 <= tFinal; t1 += dtPrint { printErr(t, y) // print intermediate result for steps := int(float64(dtPrint)/dtStep + .5); steps > 1; steps-- { y = ypStep(t, y, dtStep) t += dtStep } y = ypStep(t, y, float64(t1)-t) // adjust step to integer time t = float64(t1) } printErr(t, y) // print final result  } func printErr(t, y float64) {  fmt.Printf("y(%.1f) = %f Error: %e\n", t, y, math.Abs(actual(t)-y))  }</lang> Output: y(0.0) = 1.000000 Error: 0.000000e+00 y(1.0) = 1.562500 Error: 1.457219e-07 y(2.0) = 3.999999 Error: 9.194792e-07 y(3.0) = 10.562497 Error: 2.909562e-06 y(4.0) = 24.999994 Error: 6.234909e-06 y(5.0) = 52.562489 Error: 1.081970e-05 y(6.0) = 99.999983 Error: 1.659460e-05 y(7.0) = 175.562476 Error: 2.351773e-05 y(8.0) = 288.999968 Error: 3.156520e-05 y(9.0) = 451.562459 Error: 4.072316e-05 y(10.0) = 675.999949 Error: 5.098329e-05  ## Haskell Using GHC 7.4.1. <lang haskell>import Data.List dv :: Floating a => a -> a -> a dv = (. sqrt). (*) fy t = 1/16 * (4+t^2)^2 rk4 :: (Enum a, Fractional a)=> (a -> a -> a) -> a -> a -> a -> [(a,a)] rk4 fd y0 a h = zip ts$ scanl (flip fc) y0 ts where

 ts = [a,h ..]
fc t y = sum. (y:). zipWith (*) [1/6,1/3,1/3,1/6]
$scanl (\k f -> h * fd (t+f*h) (y+f*k)) (h * fd t y) [1/2,1/2,1]  task = mapM_ print $ map (\(x,y)-> (truncate x,y,fy x - y))
$filter (\(x,_) -> 0== mod (truncate$ 10*x) 10)
$take 101$ rk4 dv 1.0 0 0.1</lang>


Example executed in GHCi: <lang haskell>*Main> task (0,1.0,0.0) (1,1.5624998542781088,1.4572189122041834e-7) (2,3.9999990805208006,9.194792029987298e-7) (3,10.562497090437557,2.909562461184123e-6) (4,24.999993765090654,6.234909399438493e-6) (5,52.56248918030265,1.0819697635611192e-5) (6,99.99998340540378,1.6594596999652822e-5) (7,175.56247648227165,2.3517730085131916e-5) (8,288.99996843479926,3.1565204153594095e-5) (9,451.562459276841,4.0723166534917254e-5) (10,675.9999490167125,5.098330132113915e-5)</lang>

## J

Solution: <lang j>NB.*rk4 a Solve function using Runge-Kutta method NB. y is: y(ta) , ta , tb , tstep NB. u is: function to solve NB. eg: fyp rk4 1 0 10 0.1 rk4=: adverb define

'Y0 a b h'=. 4{. y
T=. a + i.@>:&.(%&h) b - a
Y=. Yt=. Y0
for_t. }: T do.
ty=. t,Yt
k1=. h * u ty
k2=. h * u ty + -: h,k1
k3=. h * u ty + -: h,k2
k4=. h * u ty + h,k3
Y=. Y, Yt=. Yt + (%6) * 1 2 2 1 +/@:* k1, k2, k3, k4
end.


T ,. Y )</lang> Example: <lang j> fy=: (%16) * [: *: 4 + *: NB. f(t,y)

  fyp=: (* %:)/                         NB. f'(t,y)
report_whole=: (10 * i. >:10)&{       NB. report at whole-numbered t values

  report_err report_whole fyp rk4 1 0 10 0.1
0       1           0
1  1.5625 _1.45722e_7
2       4 _9.19479e_7
3 10.5625 _2.90956e_6
4      25 _6.23491e_6
5 52.5625 _1.08197e_5
6     100 _1.65946e_5
7 175.562 _2.35177e_5
8     289 _3.15652e_5
9 451.562 _4.07232e_5


10 676 _5.09833e_5</lang>

Alternative solution:

The following solution replaces the for loop as well as the calculation of the increments (ks) with an accumulating suffix. <lang j>rk4=: adverb define

'Y0 a b h'=. 4{. y
T=. a + i.@>:&.(%&h) b-a
(,. [: h&(u nextY)@,/\. Y0 ,~ }.)&.|. T


)

NB. nextY a Calculate Yn+1 of a function using Runge-Kutta method NB. y is: 2-item numeric list of time t and y(t) NB. u is: function to use NB. x is: step size NB. eg: 0.001 fyp nextY 0 1 nextY=: adverb define

tableau=. 1 0.5 0.5, x * u y
ks=. (x * [: u y + (* x&,))/\. tableau
({:y) + 6 %~ +/ 1 2 2 1 * ks


)</lang>

Use:

report_err report_whole fyp rk4 1 0 10 0.1


## JavaScript

<lang JavaScript> function rk4(y, x, dx, f) {

   var k1 = dx * f(x, y),
k2 = dx * f(x + dx / 2.0,     +y + k1 / 2.0),
k3 = dx * f(x + dx / 2.0,     +y + k2 / 2.0),
k4 = dx * f(x + dx,         +y + k3);

   return y + (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0;


}

function f(t, y) {

   return t * Math.sqrt(y);


}

function actual(t) {

   return (1/16) * (t*t+4)*(t*t+4);


}

var y = 1.0,

   x = 0.0,
step = 0.1,
steps = 0,
maxSteps = 101,
sampleEveryN = 10;


while (steps < maxSteps) {

   if (steps%sampleEveryN === 0) {
console.log(x + "\t" + y + "\t" + (actual(x) - y).toExponential());
}

   y = rk4(y, x, step, f);

   // using integer math for the step addition
// to prevent floating point errors as 0.2 + 0.1 != 0.3
x = ((x * 10) + (step * 10)) / 10;
steps += 1;


} </lang>

Output:
0	1			0e+0
1	1.562499854278108	1.4572189210859676e-7
2	3.999999080520799	9.194792007782837e-7
3	10.562497090437551	2.9095624487496252e-6
4	24.999993765090636	6.234909363911356e-6
5	52.562489180302585	1.0819697415342944e-5
6	99.99998340540358	1.659459641700778e-5
7	175.56247648227125	2.3517728749311573e-5
8	288.9999684347986	3.156520142510999e-5
9	451.56245927683966	4.07231603389846e-5
10	675.9999490167097	5.098329029351589e-5


## Julia

<lang Julia> function rk4(f)

       return   (t,y,dt)->
( (dy1   )->
( (dy2   )->
( (dy3   )->
( (dy4   )->( dy1 + 2*dy2 + 2*dy3 + dy4 ) / 6
)( dt * f( t +dt  , y + dy3   ) )
)( dt * f( t +dt/2, y + dy2/2 ) )
)( dt * f( t +dt/2, y + dy1/2 ) )
)( dt * f( t      , y         ) )


end

theory(t) = (t^2 + 4.0)^2 / 16.0

tmax = 10.0 ttol = 1.e-5

t0 = 0.0 y0 = 1.0 dt = 0.1

dy = rk4( (t,y) -> t*sqrt(y) )

t = t0 y = y0

while t <= tmax

       if abs(round(t) - t) < ttol
@printf( STDOUT,"y(%4.1f)\t= %12.6f \t error: %12.6e\n",t,y,abs(y-theory(t)) )
end
y = y + dy(t,y,dt)
t = t + dt


end </lang>

Output:
y( 0.0)	=     1.000000 	 error: 0.000000e+00
y( 1.0)	=     1.562500 	 error: 1.457219e-07
y( 2.0)	=     3.999999 	 error: 9.194792e-07
y( 3.0)	=    10.562497 	 error: 2.909562e-06
y( 4.0)	=    24.999994 	 error: 6.234909e-06
y( 5.0)	=    52.562489 	 error: 1.081970e-05
y( 6.0)	=    99.999983 	 error: 1.659460e-05
y( 7.0)	=   175.562476 	 error: 2.351773e-05
y( 8.0)	=   288.999968 	 error: 3.156520e-05
y( 9.0)	=   451.562459 	 error: 4.072316e-05
y(10.0)	=   675.999949 	 error: 5.098329e-05


## Mathematica

<lang Mathematica>(* Symbolic solution *) DSolve[{y'[t] == t*Sqrt[y[t]], y[0] == 1}, y, t] Table[{t, 1/16 (4 + t^2)^2}, {t, 0, 10}]

(* Numerical solution I (not RK4) *) Table[{t, y[t], Abs[y[t] - 1/16*(4 + t^2)^2]}, {t, 0, 10}] /.

First@NDSolve[{y'[t] == t*Sqrt[y[t]], y[0] == 1}, y, {t, 0, 10}]


(* Numerical solution II (RK4) *) f[{t_, y_}] := {1, t Sqrt[y]} h = 0.1; phi[y_] := Module[{k1, k2, k3, k4},

 k1 = h*f[y];
k2 = h*f[y + 1/2 k1];
k3 = h*f[y + 1/2 k2];
k4 = h*f[y + k3];
y + k1/6 + k2/3 + k3/3 + k4/6]


solution = NestList[phi, {0, 1}, 101]; Table[{y1, y2, Abs[y2 - 1/16 (y1^2 + 4)^2]},

 {y,  solution1 ;; 101 ;; 10}]


</lang>

## MATLAB

The normally-used built-in solver is the ode45 function, which uses a non-fixed-step solver with 4th/5th order Runge-Kutta methods. The MathWorks Support Team released a package of fixed-step RK method ODE solvers on MATLABCentral. The ode4 function contained within uses a 4th-order Runge-Kutta method. Here is code that tests both ode4 and my own function, shows that they are the same, and compares them to the exact solution. <lang MATLAB>function testRK4Programs

   figure
hold on
t = 0:0.1:10;
y = 0.0625.*(t.^2+4).^2;
plot(t, y, '-k')
[tode4, yode4] = testODE4(t);
plot(tode4, yode4, '--b')
[trk4, yrk4] = testRK4(t);
plot(trk4, yrk4, ':r')
legend('Exact', 'ODE4', 'RK4')
hold off
fprintf('Time\tExactVal\tODE4Val\tODE4Error\tRK4Val\tRK4Error\n')
for k = 1:10:length(t)
fprintf('%.f\t\t%7.3f\t\t%7.3f\t%7.3g\t%7.3f\t%7.3g\n', t(k), y(k), ...
yode4(k), abs(y(k)-yode4(k)), yrk4(k), abs(y(k)-yrk4(k)))
end


end

function [t, y] = testODE4(t)

   y0 = 1;
y = ode4(@(tVal,yVal)tVal*sqrt(yVal), t, y0);


end

function [t, y] = testRK4(t)

   dydt = @(tVal,yVal)tVal*sqrt(yVal);
y = zeros(size(t));
y(1) = 1;
for k = 1:length(t)-1
dt = t(k+1)-t(k);
dy1 = dt*dydt(t(k), y(k));
dy2 = dt*dydt(t(k)+0.5*dt, y(k)+0.5*dy1);
dy3 = dt*dydt(t(k)+0.5*dt, y(k)+0.5*dy2);
dy4 = dt*dydt(t(k)+dt, y(k)+dy3);
y(k+1) = y(k)+(dy1+2*dy2+2*dy3+dy4)/6;
end


end</lang>

Output:
Time	ExactVal	ODE4Val		ODE4Error	RK4Val		RK4Error
0	  1.000		  1.000		      0		  1.000		      0
1	  1.563		  1.562		1.46e-007	  1.562		1.46e-007
2	  4.000		  4.000		9.19e-007	  4.000		9.19e-007
3	 10.563		 10.562		2.91e-006	 10.562		2.91e-006
4	 25.000		 25.000		6.23e-006	 25.000		6.23e-006
5	 52.563		 52.562		1.08e-005	 52.562		1.08e-005
6	100.000		100.000		1.66e-005	100.000		1.66e-005
7	175.563		175.562		2.35e-005	175.562		2.35e-005
8	289.000		289.000		3.16e-005	289.000		3.16e-005
9	451.563		451.562		4.07e-005	451.562		4.07e-005
10	676.000		676.000		5.1e-005	676.000		5.1e-005

## Maxima

<lang maxima>/* Here is how to solve a differential equation */ 'diff(y, x) = x * sqrt(y); ode2(%, y, x); ic1(%, x = 0, y = 1); factor(solve(%, y)); /* [y = (x^2 + 4)^2 / 16] */

/* The Runge-Kutta solver is builtin */

load(dynamics)$sol: rk(t * sqrt(y), y, 1, [t, 0, 10, 1.0])$ plot2d([discrete, sol])$/* An implementation of RK4 for one equation */ rk4(f, x0, y0, x1, n) := block([h, x, y, vx, vy, k1, k2, k3, k4],  h: bfloat((x1 - x0) / (n - 1)), x: x0, y: y0, vx: makelist(0, n + 1), vy: makelist(0, n + 1), vx[1]: x0, vy[1]: y0, for i from 1 thru n do ( k1: bfloat(h * f(x, y)), k2: bfloat(h * f(x + h / 2, y + k1 / 2)), k3: bfloat(h * f(x + h / 2, y + k2 / 2)), k4: bfloat(h * f(x + h, y + k3)), vy[i + 1]: y: y + (k1 + 2 * k2 + 2 * k3 + k4) / 6, vx[i + 1]: x: x + h ), [vx, vy]  )$

[x, y]: rk4(lambda([x, y], x * sqrt(y)), 0, 1, 10, 101)$plot2d([discrete, x, y])$

s: map(lambda([x], (x^2 + 4)^2 / 16), x)$for i from 1 step 10 thru 101 do print(x[i], " ", y[i], " ", y[i] - s[i]);</lang> ## МК-61/52 ПП 38 П1 ПП 30 П2 ПП 35 П3 2 * ПП 30 ИП2 ИП3 + 2 * + ИП1 + 3 / ИП7 + П7 П8 С/П БП 00 ИП6 ИП5 + П6 <-> ИП7 + П8 ИП8 КвКор ИП6 * ИП5 * В/О Input: 1/2 (h/2) - Р5, 1 (y0) - Р8 and Р7, 0 (t0) - Р6. ## OCaml <lang ocaml>let y' t y = t *. sqrt y let exact t = let u = 0.25*.t*.t +. 1.0 in u*.u let rk4_step (y,t) h =  let k1 = h *. y' t y in let k2 = h *. y' (t +. 0.5*.h) (y +. 0.5*.k1) in let k3 = h *. y' (t +. 0.5*.h) (y +. 0.5*.k2) in let k4 = h *. y' (t +. h) (y +. k3) in (y +. (k1+.k4)/.6.0 +. (k2+.k3)/.3.0, t +. h)  let rec loop h n (y,t) =  if n mod 10 = 1 then Printf.printf "t = %f,\ty = %f,\terr = %g\n" t y (abs_float (y -. exact t)); if n < 102 then loop h (n+1) (rk4_step (y,t) h)  let _ = loop 0.1 1 (1.0, 0.0)</lang> Output: t = 0.000000, y = 1.000000, err = 0 t = 1.000000, y = 1.562500, err = 1.45722e-07 t = 2.000000, y = 3.999999, err = 9.19479e-07 t = 3.000000, y = 10.562497, err = 2.90956e-06 t = 4.000000, y = 24.999994, err = 6.23491e-06 t = 5.000000, y = 52.562489, err = 1.08197e-05 t = 6.000000, y = 99.999983, err = 1.65946e-05 t = 7.000000, y = 175.562476, err = 2.35177e-05 t = 8.000000, y = 288.999968, err = 3.15652e-05 t = 9.000000, y = 451.562459, err = 4.07232e-05 t = 10.000000, y = 675.999949, err = 5.09833e-05 ## Octave  This example is incomplete. Please ensure that it meets all task requirements and remove this message. Implementing Runge-Kutta in octave is a bit of a hassle. For now we'll showcase how to solve an ordinary differential equation using the lsode function. <lang octave>function ydot = f(y, t)  ydot = t * sqrt( y );  endfunction t = [0:10]'; y = lsode("f", 1, t); [ t, y, y - 1/16 * (t.**2 + 4).**2 ]</lang> Output: ans = 0.00000 1.00000 0.00000 1.00000 1.56250 0.00000 2.00000 4.00000 0.00000 3.00000 10.56250 0.00000 4.00000 25.00000 0.00000 5.00000 52.56250 0.00000 6.00000 100.00000 0.00000 7.00000 175.56250 0.00000 8.00000 289.00001 0.00001 9.00000 451.56251 0.00001 10.00000 676.00001 0.00001 ## PARI/GP Translation of: C <lang parigp>rk4(f,dx,x,y)={  my(k1=dx*f(x,y), k2=dx*f(x+dx/2,y+k1/2), k3=dx*f(x+dx/2,y+k2/2), k4=dx*f(x+dx,y+k3)); y + (k1 + 2*k2 + 2*k3 + k4) / 6  }; rate(x,y)=x*sqrt(y); go()={  my(x0=0,x1=10,dx=.1,n=1+(x1-x0)\dx,y=vector(n)); y[1]=1; for(i=2,n,y[i]=rk4(rate, dx, x0 + dx * (i - 1), y[i-1])); print("x\ty\trel. err.\n------------"); forstep(i=1,n,10, my(x=x0+dx*i,y2=(x^2/4+1)^2); print(x "\t" y[i] "\t" y[i]/y2 - 1) )  }; go()</lang> Output: x y rel. err. ------------ 0.100000000 1 -0.00498131231 1.10000000 1.68999982 -0.00383519474 2.10000000 4.40999894 -0.00237694942 3.10000000 11.5599968 -0.00146924588 4.10000000 27.0399933 -0.000961094862 5.10000000 56.2499884 -0.000666538719 6.10000000 106.089982 -0.000485427212 7.10000000 184.959975 -0.000367681962 8.10000000 302.759966 -0.000287408941 9.10000000 470.889955 -0.000230470905 ## Pascal Translation of: Ada This code has been compiled using Free Pascal 2.6.2. <lang pascal>program RungeKuttaExample; uses sysutils; type  TDerivative = function (t, y : Real) : Real;  procedure RungeKutta(yDer : TDerivative;  var t, y : array of Real; dt : Real);  var  dy1, dy2, dy3, dy4 : Real; idx : Cardinal;  begin  for idx := Low(t) to High(t) - 1 do begin dy1 := dt * yDer(t[idx], y[idx]); dy2 := dt * yDer(t[idx] + dt / 2.0, y[idx] + dy1 / 2.0); dy3 := dt * yDer(t[idx] + dt / 2.0, y[idx] + dy2 / 2.0); dy4 := dt * yDer(t[idx] + dt, y[idx] + dy3); t[idx + 1] := t[idx] + dt; y[idx + 1] := y[idx] + (dy1 + 2.0 * (dy2 + dy3) + dy4) / 6.0; end;  end; function CalcError(t, y : Real) : Real; var  trueVal : Real;  begin  trueVal := sqr(sqr(t) + 4.0) / 16.0; CalcError := abs(trueVal - y);  end; procedure Print(t, y : array of Real;  modnum : Integer);  var  idx : Cardinal;  begin  for idx := Low(t) to High(t) do begin if idx mod modnum = 0 then begin WriteLn(Format('y(%4.1f) = %12.8f Error: %12.6e', [t[idx], y[idx], CalcError(t[idx], y[idx])])); end; end;  end; function YPrime(t, y : Real) : Real; begin  YPrime := t * sqrt(y);  end; const  dt = 0.10; N = 100;  var  tArr, yArr : array [0..N] of Real;  begin  tArr[0] := 0.0; yArr[0] := 1.0; RungeKutta(@YPrime, tArr, yArr, dt); Print(tArr, yArr, 10);  end.</lang> Output: y( 0.0) = 1.00000000 Error: 0.00000E+000 y( 1.0) = 1.56249985 Error: 1.45722E-007 y( 2.0) = 3.99999908 Error: 9.19479E-007 y( 3.0) = 10.56249709 Error: 2.90956E-006 y( 4.0) = 24.99999377 Error: 6.23491E-006 y( 5.0) = 52.56248918 Error: 1.08197E-005 y( 6.0) = 99.99998341 Error: 1.65946E-005 y( 7.0) = 175.56247648 Error: 2.35177E-005 y( 8.0) = 288.99996843 Error: 3.15652E-005 y( 9.0) = 451.56245928 Error: 4.07232E-005 y(10.0) = 675.99994902 Error: 5.09833E-005  ## Perl There are many ways of doing this. Here we define the runge_kutta function as a function of ${\displaystyle y'}$ and ${\displaystyle \delta t}$, returning a closure which itself takes ${\displaystyle (t,y)}$ as argument and returns the next ${\displaystyle (t,y)}$. Notice how we have to use sprintf to deal with floating point rounding. See perlfaq4. <lang perl>sub runge_kutta {  my ($yp, $dt) = @_; sub {  my ($t, $y) = @_; my @dy =$dt * $yp->($t , $y ); push @dy,$dt * $yp->($t + $dt/2,$y + $dy[0]/2 ); push @dy,$dt * $yp->($t + $dt/2,$y + $dy[1]/2 ); push @dy,$dt * $yp->($t + $dt ,$y + $dy[2] ); return$t + $dt,$y + ($dy[0] + 2*$dy[1] + 2*$dy[2] +$dy[3]) / 6;

   }


}

my $RK = runge_kutta sub {$_[0] * sqrt $_[1] }, .1; for(  my ($t, $y) = (0, 1); sprintf("%.0f",$t) <= 10;
($t,$y) = $RK->($t, $y)  ) {  printf "y(%2.0f) = %12f ± %e\n",$t, $y, abs($y - ($t**2 + 4)**2 / 16) if sprintf("%.4f",$t) =~ /0000$/;  }</lang> 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 ## Perl 6 <lang perl6>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 + δt, $y + δy($t, $y, δt))  ) {  printf "y(%2d) = %12f ± %e\n",$t, $y, abs($y - ($t**2 + 4)**2 / 16) if$t.narrow ~~ Int;


}</lang>

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

## PL/I

<lang PL/I> Runge_Kutta: procedure options (main); /* 10 March 2014 */

  declare (y, dy1, dy2, dy3, dy4) float (18);
declare t fixed decimal (10,1);
declare dt float (18) static initial (0.1);

  y = 1;
do t = 0 to 10 by 0.1;
dy1 = dt * ydash(t, y);
dy2 = dt * ydash(t + dt/2, y + dy1/2);
dy3 = dt * ydash(t + dt/2, y + dy2/2);
dy4 = dt * ydash(t + dt,   y + dy3);

     if mod(t, 1.0) = 0 then
put skip edit('y(', trim(t), ')=', y, ', error = ', abs(y - (t**2 + 4)**2 / 16 ))
(3 a, column(9), f(16,10), a, f(13,10));
y = y + (dy1 + 2*dy2 + 2*dy3 + dy4)/6;
end;


ydash: procedure (t, y) returns (float(18));

  declare (t, y) float (18) nonassignable;
return ( t*sqrt(y) );


end ydash;

end Runge_kutta; </lang> Output:-

y(0.0)=     1.0000000000, error =  0.0000000000
y(1.0)=     1.5624998543, error =  0.0000001457
y(2.0)=     3.9999990805, error =  0.0000009195
y(3.0)=    10.5624970904, error =  0.0000029096
y(4.0)=    24.9999937651, error =  0.0000062349
y(5.0)=    52.5624891803, error =  0.0000108197
y(6.0)=    99.9999834054, error =  0.0000165946
y(7.0)=   175.5624764823, error =  0.0000235177
y(8.0)=   288.9999684348, error =  0.0000315652
y(9.0)=   451.5624592768, error =  0.0000407232
y(10.0)=  675.9999490167, error =  0.0000509833


## Python

<lang Python>def RK4(f):

   return lambda t, y, dt: (
lambda dy1: (
lambda dy2: (
lambda dy3: (
lambda dy4: (dy1 + 2*dy2 + 2*dy3 + dy4)/6
)( dt * f( t + dt  , y + dy3   ) )


)( dt * f( t + dt/2, y + dy2/2 ) ) )( dt * f( t + dt/2, y + dy1/2 ) ) )( dt * f( t , y ) )

def theory(t): return (t**2 + 4)**2 /16

from math import sqrt dy = RK4(lambda t, y: t*sqrt(y))

t, y, dt = 0., 1., .1 while t <= 10:

   if abs(round(t) - t) < 1e-5:


print("y(%2.1f)\t= %4.6f \t error: %4.6g" % ( t, y, abs(y - theory(t))))

   t, y = t + dt, y + dy( t, y, dt )


</lang>

Output:
y(0.0)	= 1.000000 	 error:    0
y(1.0)	= 1.562500 	 error: 1.45722e-07
y(2.0)	= 3.999999 	 error: 9.19479e-07
y(3.0)	= 10.562497 	 error: 2.90956e-06
y(4.0)	= 24.999994 	 error: 6.23491e-06
y(5.0)	= 52.562489 	 error: 1.08197e-05
y(6.0)	= 99.999983 	 error: 1.65946e-05
y(7.0)	= 175.562476 	 error: 2.35177e-05
y(8.0)	= 288.999968 	 error: 3.15652e-05
y(9.0)	= 451.562459 	 error: 4.07232e-05
y(10.0)	= 675.999949 	 error: 5.09833e-05

## Ruby

<lang ruby> def calc_rk4(f)

 return ->(t,y,dt){
->(dy1   ){
->(dy2   ){
->(dy3   ){
->(dy4   ){ ( dy1 + 2*dy2 + 2*dy3 + dy4 ) / 6 }.call(
dt * f.call( t + dt  , y + dy3   ))}.call(
dt * f.call( t + dt/2, y + dy2/2 ))}.call(
dt * f.call( t + dt/2, y + dy1/2 ))}.call(
dt * f.call( t       , y         ))}


end

TIME_MAXIMUM, WHOLE_TOLERANCE = 10.0, 1.0e-5 T_START, Y_START, DT = 0.0, 1.0, 0.10

def my_diff_eqn(t,y) ; t * Math.sqrt(y)  ; end def my_solution(t ) ; (t**2 + 4)**2 / 16  ; end def find_error(t,y) ; (y - my_solution(t)).abs  ; end def is_whole?(t ) ; (t.round - t).abs < WHOLE_TOLERANCE ; end

dy = calc_rk4( ->(t,y){my_diff_eqn(t,y)} )

t, y = T_START, Y_START while t <= TIME_MAXIMUM

 printf("y(%4.1f)\t= %12.6f \t error: %12.6e\n",t,y,find_error(t,y)) if is_whole?(t)
t, y = t + DT, y + dy.call(t,y,DT)


end </lang>

Output:
y( 0.0)	=     1.000000 	 error: 0.000000e+00
y( 1.0)	=     1.562500 	 error: 1.457219e-07
y( 2.0)	=     3.999999 	 error: 9.194792e-07
y( 3.0)	=    10.562497 	 error: 2.909562e-06
y( 4.0)	=    24.999994 	 error: 6.234909e-06
y( 5.0)	=    52.562489 	 error: 1.081970e-05
y( 6.0)	=    99.999983 	 error: 1.659460e-05
y( 7.0)	=   175.562476 	 error: 2.351773e-05
y( 8.0)	=   288.999968 	 error: 3.156520e-05
y( 9.0)	=   451.562459 	 error: 4.072316e-05
y(10.0)	=   675.999949 	 error: 5.098329e-05


## Racket

See Euler method#Racket for implementation of simple general ODE-solver.

The Runge-Kutta method <lang racket> (define (RK4 F δt)

 (λ (t y)
(define δy1 (* δt (F t y)))
(define δy2 (* δt (F (+ t (* 1/2 δt)) (+ y (* 1/2 δy1)))))
(define δy3 (* δt (F (+ t (* 1/2 δt)) (+ y (* 1/2 δy2)))))
(define δy4 (* δt (F (+ t δt) (+ y δy1))))
(list (+ t δt)
(+ y (* 1/6 (+ δy1 (* 2 δy2) (* 2 δy3) δy4))))))


</lang>

The method modifier which divides each time-step into n sub-steps: <lang racket> (define ((step-subdivision n method) F h)

 (λ (x . y) (last (ODE-solve F (cons x y)
#:x-max (+ x h)
#:step (/ h n)
#:method method))))


</lang>

Usage: <lang racket> (define (F t y) (* t (sqrt y)))

(define (exact-solution t) (* 1/16 (sqr (+ 4 (sqr t)))))

(define numeric-solution

   (ODE-solve F '(0 1) #:x-max 10 #:step 1 #:method (step-subdivision 10 RK4)))


(for ([s numeric-solution])

 (match-define (list t y) s)
(printf "t=~a\ty=~a\terror=~a\n" t y (- y (exact-solution t))))


</lang>

Output:
t=0	y=1	                error=0
t=1	y=1.562499854278108	error=-1.4572189210859676e-07
t=2	y=3.999999080520799	error=-9.194792007782837e-07
t=3	y=10.562497090437551	error=-2.9095624487496252e-06
t=4	y=24.999993765090636	error=-6.234909363911356e-06
t=5	y=52.562489180302585	error=-1.0819697415342944e-05
t=6	y=99.99998340540358	error=-1.659459641700778e-05
t=7	y=175.56247648227125	error=-2.3517728749311573e-05
t=8	y=288.9999684347986	error=-3.156520142510999e-05
t=9	y=451.56245927683966	error=-4.07231603389846e-05
t=10	y=675.9999490167097	error=-5.098329029351589e-05


Graphical representation:

<lang racket> > (require plot) > (plot (list (function exact-solution 0 10 #:label "Exact solution")

             (points numeric-solution #:label "Runge-Kutta method"))
#:x-label "t" #:y-label "y(t)")


</lang>

## REXX

<lang rexx>/*REXX program uses the Runge-Kutta method to solve the differential */ /* ____ */ /*equation: y'(t)=t²√y(t) which has the exact solution: y(t)=(t²+4)²/16*/

numeric digits 40; d=digits()%2 /*use forty digits, show ½ that. */ x0=0; x1=10; dx=.1; n=1 + (x1-x0) / dx; y.=1

      do m=1  for n-1;                mm=m-1
y.m=Runge_Kutta(dx, x0+dx*mm, y.mm)
end   /*m*/


say center(x,13,'─') center(y,d,'─') ' ' center('relative error',d,'─')

      do i=0  to n-1  by 10;         x=(x0+dx*i)/1;       y2=(x*x/4+1)**2
relE=format(y.i/y2-1,,13)/1;   if relE=0  then relE=' 0'
say  center(x,13)   right(format(y.i,,12),d)    '  '   left(relE,d)
end   /*i*/


exit /*stick a fork in it, we're done.*/ /*──────────────────────────────────RATE subroutine─────────────────────*/ rate: return arg(1)*sqrt(arg(2)) /*──────────────────────────────────Runge_Kutta subroutine──────────────*/ Runge_Kutta: procedure; parse arg dx,x,y

                                         k1 = dx * rate(x      , y      )
k2 = dx * rate(x+dx/2 , y+k1/2 )
k3 = dx * rate(x+dx/2 , y+k2/2 )
k4 = dx * rate(x+dx   , y+k3   )


return y + (k1 + 2*k2 + 2*k3 + k4) / 6 /*──────────────────────────────────SQRT subroutine─────────────────────*/ sqrt: procedure; parse arg x; if x=0 then return 0; d=digits()

      numeric digits 11;        g=.sqrtG()
do j=0 while p>9; m.j=p; p=p%2+1; end;  do k=j+5 to 0 by -1
if m.k>11  then numeric digits m.k
g=.5*(g+x/g); end;        numeric digits d;     return g/1


.sqrtG: numeric form; m.=11; p=d+d%4+2

      parse value format(x,2,1,,0) 'E0' with g 'E' _ .; return g*.5'E'_%2</lang>


output

──────X────── ─────────Y──────────   ───relative error───
0             1.000000000000     0
1             1.562499854278    -9.3262010934636E-8
2             3.999999080521    -2.2986980018794E-7
3            10.562497090438    -2.7546153355698E-7
4            24.999993765091    -2.4939637458826E-7
5            52.562489180303    -2.058444217357E-7
6            99.999983405404    -1.6594596403363E-7
7           175.562476482271    -1.3395644712797E-7
8           288.999968434799    -1.092221504E-7
9           451.562459276840    -9.018277747572E-8
10           675.999949016709    -7.5419068845528E-8


## Run BASIC

<lang Runbasic>y = 1 while t <= 10

  k1	=  t        * sqr(y)
k2	= (t + .05) * sqr(y + .05 * k1)
k3	= (t + .05) * sqr(y + .05 * k2)
k4	= (t + .1)  * sqr(y + .1  * k3)


if right$(using("##.#",t),1) = "0" then print "y(";using("##",t);") ="; using("####.#######", y);chr$(9);"Error ="; (((t^2 + 4)^2) /16) -y

   y = y + .1 *(k1 + 2 * (k2 + k3) + k4) / 6
t = t + .1


wend end</lang>

Output:
y( 0) =   1.0000000	Error =0
y( 1) =   1.5624999	Error =1.45721892e-7
y( 2) =   3.9999991	Error =9.19479203e-7
y( 3) =  10.5624971	Error =2.90956246e-6
y( 4) =  24.9999938	Error =6.23490939e-6
y( 5) =  52.5624892	Error =1.08196973e-5
y( 6) =  99.9999834	Error =1.65945961e-5
y( 7) = 175.5624765	Error =2.3517728e-5
y( 8) = 288.9999684	Error =3.15652e-5
y( 9) = 451.5624593	Error =4.07231581e-5
y(10) = 675.9999490	Error =5.09832864e-5


## Standard ML

<lang sml>fun step y' (tn,yn) dt =

   let
val dy1 = dt * y'(tn,yn)
val dy2 = dt * y'(tn + 0.5 * dt, yn + 0.5 * dy1)
val dy3 = dt * y'(tn + 0.5 * dt, yn + 0.5 * dy2)
val dy4 = dt * y'(tn + dt, yn + dy3)
in
(tn + dt, yn + (1.0 / 6.0) * (dy1 + 2.0*dy2 + 2.0*dy3 + dy4))
end


(* Suggested test case *) fun testy' (t,y) =

   t * Math.sqrt y


fun testy t =

   (1.0 / 16.0) * Math.pow(Math.pow(t,2.0) + 4.0, 2.0)


(* Test-runner that iterates the step function and prints the results. *) fun test t0 y0 dt steps print_freq y y' =

   let
fun loop i (tn,yn) =
if i = steps then ()
else
let
val (t1,y1) = step y' (tn,yn) dt
val y1' = y tn
val () = if i mod print_freq = 0 then
(print ("Time: " ^ Real.toString tn ^ "\n");
print ("Exact: " ^ Real.toString y1' ^ "\n");
print ("Approx: " ^ Real.toString yn ^ "\n");
print ("Error: " ^ Real.toString (y1' - yn) ^ "\n\n"))
else ()
in
loop (i+1) (t1,y1)
end
in
loop 0 (t0,y0)
end


(* Run the suggested test case *) val () = test 0.0 1.0 0.1 101 10 testy testy'</lang> Output

Time: 0.0
Exact: 1.0
Approx: 1.0
Error: ~1.11022302463E~16

Time: 1.0
Exact: 1.5625
Approx: 1.56249985428
Error: 1.45722452549E~07

Time: 2.0
Exact: 4.0
Approx: 3.99999908052
Error: 9.19479203443E~07

Time: 3.0
Exact: 10.5625
Approx: 10.5624970904
Error: 2.90956245586E~06

Time: 4.0
Exact: 25.0
Approx: 24.9999937651
Error: 6.23490938878E~06

Time: 5.0
Exact: 52.5625
Approx: 52.5624891803
Error: 1.08196973727E~05

Time: 6.0
Exact: 100.0
Approx: 99.9999834054
Error: 1.65945961186E~05

Time: 7.0
Exact: 175.5625
Approx: 175.562476482
Error: 2.35177280956E~05

Time: 8.0
Exact: 289.0
Approx: 288.999968435
Error: 3.15651997767E~05

Time: 9.0
Exact: 451.5625
Approx: 451.562459277
Error: 4.07231581221E~05

Time: 10.0
Exact: 676.0
Approx: 675.999949017
Error: 5.09832866555E~05

## Tcl

<lang tcl>package require Tcl 8.5

1. Hack to bring argument function into expression

proc tcl::mathfunc::dy {t y} {upvar 1 dyFn dyFn; $dyFn$t $y} proc rk4step {dyFn y* t* dt} {  upvar 1${y*} y ${t*} t set dy1 [expr {$dt * dy($t,$y)}]
set dy2 [expr {$dt * dy($t+$dt/2,$y+$dy1/2)}] set dy3 [expr {$dt * dy($t+$dt/2, $y+$dy2/2)}]
set dy4 [expr {$dt * dy($t+$dt,$y+$dy3)}] set y [expr {$y + ($dy1 + 2*$dy2 + 2*$dy3 +$dy4)/6.0}]
set t [expr {$t +$dt}]


}

proc y {t} {expr {($t**2 + 4)**2 / 16}} proc δy {t y} {expr {$t * sqrt($y)}} proc printvals {t y} {  set err [expr {abs($y - [y $t])}] puts [format "y(%.1f) = %.8f\tError: %.8e"$t $y$err]


}

set t 0.0 set y 1.0 set dt 0.1 printvals $t$y for {set i 1} {$i <= 101} {incr i} {  rk4step δy y t$dt
if {$i%10 == 0} {  printvals$t \$y

   }


}</lang>

Output:
y(0.0) = 1.00000000	Error: 0.00000000e+00
y(1.0) = 1.56249985	Error: 1.45721892e-07
y(2.0) = 3.99999908	Error: 9.19479203e-07
y(3.0) = 10.56249709	Error: 2.90956245e-06
y(4.0) = 24.99999377	Error: 6.23490939e-06
y(5.0) = 52.56248918	Error: 1.08196973e-05
y(6.0) = 99.99998341	Error: 1.65945961e-05
y(7.0) = 175.56247648	Error: 2.35177280e-05
y(8.0) = 288.99996843	Error: 3.15652000e-05
y(9.0) = 451.56245928	Error: 4.07231581e-05
y(10.0) = 675.99994902	Error: 5.09832864e-05

## zkl

Translation of: OCaml

<lang zkl>fcn yp(t,y) { t * y.sqrt() } fcn exact(t){ u:=0.25*t*t + 1.0; u*u }

fcn rk4_step([(y,t)],h){

  k1:=h * yp(t,y);
k2:=h * yp(t + 0.5*h, y + 0.5*k1);
k3:=h * yp(t + 0.5*h, y + 0.5*k2);
k4:=h * yp(t + h, y + k3);
T(y + (k1+k4)/6.0 + (k2+k3)/3.0, t + h);


}

fcn loop(h,n,[(y,t)]){

  if(n % 10 == 1)
print("t = %f,\ty = %f,\terr = %g\n".fmt(t,y,(y - exact(t)).abs()));
if(n < 102) return(loop(h,(n+1),rk4_step(T(y,t),h))) //tail recursion


}</lang>

Output:
loop(0.1,1,T(1.0, 0.0))
t = 0.000000,	y = 1.000000,	err = 0
t = 1.000000,	y = 1.562500,	err = 1.45722e-07
t = 2.000000,	y = 3.999999,	err = 9.19479e-07
t = 3.000000,	y = 10.562497,	err = 2.90956e-06
t = 4.000000,	y = 24.999994,	err = 6.23491e-06
t = 5.000000,	y = 52.562489,	err = 1.08197e-05
t = 6.000000,	y = 99.999983,	err = 1.65946e-05
t = 7.000000,	y = 175.562476,	err = 2.35177e-05
t = 8.000000,	y = 288.999968,	err = 3.15652e-05
t = 9.000000,	y = 451.562459,	err = 4.07232e-05
t = 10.000000,	y = 675.999949,	err = 5.09833e-05