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 }




 
 [Demo of EDSAC library subroutine G1: Runge-Kutta solution of differential equations.
  Full description is in Wilkes, Wheeler & Gill, 1951 edn, pages 32-34, 86-87, 132-134.
 
  Before using G1, we need to fix n, m, a, b, c, d, as defined in WWG pages 86-87:
  n = number of equations (2 for the Rosetta Code example).
  2^m = multiplier for the hy', as large as possible without causing numeric overflow;
    with the scaling chosen here, m = 5.
  Variables y are stored in n consecutive long locations, the last of which is aD.
  Scaled derivatives (2^m)hy' in n consecutive long locations, the last of which is bD.
  G1 uses working variables in n consecutive long locations, the last of which is cD.
  d = address of user-supplied auxiliary subroutine, which calculates the (2^m)hy'.
 
  For convenience, keep G1 and its storage together. Start at (say) 400 and place:
      variables y at 400D, 402D;
      scaled derivatives at 404D, 406D;
      workspace for G1 at 408D, 410D;
      G1 itself at 412.
  If the base address is placed in location 51 at load time, all the above
    addresses can be accessed via the G parameter:]
            T   51 K
            P  400 F
 [Now set up the 6 preset parameters specified in WWG:]
            T   45 K
            P    2#G [H parameter: P a D]
            P    4 F [N parameter: P 2n F]
            P    4 F [M parameter: P (b-a) F, or V (2048-a+b) F if a > b]
            P    4 F [& parameter: P (c-b) F, or V (2048-b+c) F if b > c]
            P    8 F [L parameter: P 2^(m-2) F]
            P  300 F [X parameter: P d F]
 [For other addresses in the program we can optionally use some more parameters:]
            T   52 K
            P  120 F [A parameter: main routine]
            P   56 F [B parameter: print subroutine P1 from EDSAC library]
            P  350 F [C parameter: constants for Rosetta code example]
            P   78 F [V parameter: square root subroutine]
 
 [Library subroutine to read constants; runs at load time and is then overwritten.
  R5, for decimal fractions, seems to be unavailable (lost?), so the values are
    here read in as 35-bit integers (i.e. times 2^34) by R2.
  Values are: 0.001, initial value of y
              (2^23)/(10^7) and 25/(2^10) for use in calculations
              0.5/(10^9) for rounding to 9 d.p. (print routine P1 doesn't do this)]
    GKT20FVDL8FA40DUDTFI40FA40FS39FG@S2FG23FA5@T5@E4@E13Z
            T#C
    17179869F14411518808F419430400F9#
            TZ
 
 [Library subroutine M3; prints header at load time and is then overwritten.]
    PFGKIFAFRDLFUFOFE@A6FG@E8FEZPF
    *SCALED!FOR!EDSAC@&!!TIME!!!!!!!!!Y!VIA!RK!!!!!Y!DIRECT@&
            ....PK  [end text with some blank tape]
 
 [Runge-Kutta: auxiliary subroutine to calculate (2^m)*h*(dy1/dt) and (2^m)*h*(dy2/dt)
    from y1, y2, where y1 is the function y in Rosetta Code (but scaled) and y2 = t.
  For the Rosetta code example we're using  m = 5, h = 2^(-7)]
          E25K  TX   GK
          A3F  T20@           [set up return as usual]
          H2#G  V2#G  TD      [acc := t^2, temp store in 0D]
          H#G  VD  LD  YF  TD [y1 times t^2, shift left, round, temp store in 0D]
          H2#C  VD  YF  T4D   [times (2^23)/(10^7), round, to 4D for square root]
     [14] A14@  GV  A4D  T4#G [call square root, result in 4D, copy to (2^m)hy']
          A21@  T6#G          [1/4, i.e. (2^m)h with m and h as above, to (2^m)ht']
     [20] ZF                  [overwritten by jump back to caller]
     [21] RF                  [constant 1/4]
 
  [Main routine, with two subroutines in the same address block as the main routine.]
          E25K  TA  GK
      [0] #F                  [figures shift on teleprinter]
      [1] MF                  [decimal point (in figures mode)]
      [2] !F  @F  &F          [space, carriage return, line feed,]
      [5] K4096F              [null char]
      [6] P100F               [constant: nr of Runge-Kutta steps (in address field)]
      [7] PF                  [negative count of Runge-Kutta steps]
      [8] P10F                [constant: number of steps between printed values]
      [9] PF                  [negative count of steps between printed values]
     [Enter with acc = 0]
     [10] O@                  [set teleprinter to figures]
          S6@  T7@            [init negative count of R-K steps]
          S8@  T9@            [init negative count of print steps]
     [Before using library subroutine G1, clear its working registers (WWG page 33)]
          T8#G  T10#G
     [Set up initial values of y1 and y2 (where y2 = t)]
          A#C  T#G            [load 0.001 from constants section, store in y1]
          T2#G                [y2 = t = 0]
     [20] A20@  G40@          [call subroutine to print initial values]
     [Loop round Runge-Kutta steps]
     [22] TF  A23@  G12G      [clear accumulator, call G1 for Runge-Kutta step]
          A9@  A2F  U9@       [update negative print count]
          G33@                [skip printing if not reached 0]
          S8@  T9@            [reset negative print count]
          A31@  G40@          [call subroutine to print values]
     [33] TF                  [clear accumulator]
          A7@  A2F  U7@       [increment negative count of Runge-Kutta steps]
          G22@                [loop till count = 0]
          O5@  ZF             [flush teleprinter buffer; stop]
 
 [Subroutine to print y1 as calculated (1) by Runge-Kutta (2) direct from formula]
     [40] A3F  T71@           [set up return as usual]
          A2#G  TD            [latest t (= y2) from Runge-Kutta, to 0D for printing]
     [44] A44@  G72@          [call subroutine to print t]
          O2@  O2@            [followed by 2 spaces]
          A#G  TD             [latest y1 from Runge-Kutta, to 0D for printing]
     [50] A50@  G72@          [call subroutine to print y1]
          O2@  O2@            [followed by 2 spaces]
          A    4#C            [load constant 25/(2^10)]
          H2#G  V2#G  TD      [add t^2, temp store result in 0D]
          HD  VD  LD  YF  TD  [square, shift 1 left, round, result to 0D]
          H2#C  VD  YF  TD    [times (2^23)/(10^7), round, to 0D for printing]
     [67] A67@  G72@          [call subroutine to print y]
          O3@  O4@            [print CR, LF]
     [71] ZF                  [overwritten by jump back to caller]
 
 [Second-level subroutine to print number in 0D to 9 decimal places]
     [72] A3F  T82@           [set up return as usual]
          AD  A6#C  TD        [load number, add decimal rounding, to 0D for printing]
          O81@  O1@           [print '0.' since P1 doesn't do so]
          A79@  GB            [call library subroutine P1 for printing]
     [81] P9F                 [parameter for P1, 9 decimals]
     [82] ZF                  [overwritten by jump back to caller]
 
 [Library subroutine G1 for Runge-Kutta process. 66 locations, even address.]
            E25K  T12G
    GKT4#ZH682DT6#ZPNT12#Z!1405DT14#ZTHT16#ZT2HTZA3FT61@A31@G63@&FT6ZPN
    T8ZMMO&H4@A20@E23@T14ZAHT16ZA2HT18ZH12#@S12#@T12#@E28@H4#@T4DUFS38@
    A25@T38@S6#@A16#@U46#@A8@U37@A9@U55@A24@T39@ZFR1057#@ZFYFU6DV6DRLYF
    UDZFZFADLDADLLS6DN4DYFZFA46#@S14#@G29@A65@S11@ZFA35@U65@GXZF
 
 [Replacement for library routine S2 (square root). 38 locations, even address.
  Advantages: More accurate for small values of the argument.
              Calculates sqrt(0) without going into an infinite loop.
  Disadvantages: Longer and slower than S2 (calculates one bit at a time).]
            E25K  TV
    GKA3FT31@A4DG32@A33@T36#@T4DA33@RDU34#@RDS4DS33@A36#@G22@T36#@A4DS34#@
    T4DA36#@A33@G25@TFA36#@S33@A36#@T36#@A34#@RDYFG9@ZFZFK4096FPFPFPFPF
 
 [Library subroutine P1 - print a single positive number. 21 locations.
  Prints number in 0D to n places of decimals, where
  n is specified by 'P n F' pseudo-order after subroutine call.]
            E25K  TB
    GKA18@U17@S20@T5@H19@PFT5@VDUFOFFFSFL4FTDA5@A2FG6@EFU3FJFM1F
 
  [Define entry point in main routine]
            E25K  TA  GK
            E10Z  PF  [enter at relative address 10 with accumulator = 0]
 