1.

program flip
! test coin flipping
implicit none
  integer:: i, nrofcoins, heads=0
  real :: rand

  write(*,'("nr of coins:")', advance='no')
  read (*,*) nrofcoins

  do i=1,nrofcoins
    call random_number(rand)
    if (rand < 0.5) heads = heads+1
  end do

  write(*,'("frequency of heads:", F8.5)') (1.0*heads)/nrofcoins

end program 

> ./a.out
nr of coins:10
frequency of heads: 0.60000

> ./a.out
nr of coins:100
frequency of heads: 0.42000

> ./a.out
nr of coins:1000
frequency of heads: 0.48500

> ./a.out
nr of coins:10000
frequency of heads: 0.49240

> ./a.out
nr of coins:1000000 
frequency of heads: 0.50001


2.
program euler
! test the explicit Euler method by solving y'=x-y
  implicit none
  real, dimension(4) :: h=(/0.1, 0.01, 0.001, 0.0001 /)
  integer i
  ! try different step sizes
  do i=1,4
    call euler2(0.0, 1.0, 0.0, h(i))
  end do
contains

real function f(x, y)
  real, intent(in) :: x, y
  f = x-y
end function

subroutine euler1 (a, b, y0, h)
! explicit Euler
! solve the differential equation y'=f(x,y)
! with the initial value y(a)=y0 in the range a <= x <= b
  real, intent(in) :: a, b, y0, h
  real :: x, y, dy,  e
  x=a
  y=y0
  dy=f(x,y)
  do while (x < b-0.5*h)
    x=x+h
    y=y+h*dy
    dy=f(x,y)
  end do

  e=exp(-b)+b-1
  e=abs((e-y)/e)
  write(*,'(4F10.6)') h,x,y,e
end subroutine 
end program

  0.100000  1.000000  0.348678  0.052194
  0.010000  0.999999  0.366032  0.005021
  0.001000  0.999991  0.367693  0.000507
  0.000100  0.999954  0.367802  0.000210



3.
program euler
! test the implicit Euler method by solving y'=x-y
  implicit none
  real, dimension(4) :: h=(/0.1, 0.01, 0.001, 0.0001 /)
  integer i
  ! try different step sizes
  do i=1,4
    call euler2(0.0, 1.0, 0.0, h(i))
  end do
contains

real function f(x, y)
  real, intent(in) :: x, y
  f = x-y
end function

subroutine euler2 (a, b, y0, h)
! implicit Euler
! solve the differential equation y'=f(x,y)
! with the initial value y(a)=y0 in the range a <= x <= b
  real, intent(in) :: a, b, y0, h
  real :: x, y, y1, dy, dy0, e
  x=a
  y=y0
  y1=y
  dy0=f(x,y)
  do while (x < b-h+0.1*h)
    x=x+h
    ! predictor
    y=y1+h*dy0/2
    dy=f(x,y)
    ! corrector
    y=y1+h*(dy0+dy)/2
    dy0=dy
    y1=y
  end do

  e=exp(-b)+b-1
  e=abs((e-y)/e)
  write(*,'(4F10.6)') h,x,y,e
end subroutine 
end program

  0.100000  1.000000  0.381245  0.036330
  0.010000  0.999999  0.369202  0.003595
  0.001000  0.999991  0.368009  0.000353
  0.000100  0.999954  0.367834  0.000123


4.
program runge
! test the RK4 by solving y'=x-y
  implicit none
  real, dimension(4) :: h=(/0.1, 0.01, 0.001, 0.0001 /)
  integer i
  do i=1,4
    call rk(0.0, 1.0, 0.0, h(i))
  end do
contains
real function f(x,y)
  real, intent(in) :: x,y
  f=x-y
end function
subroutine rk(a, b, y0, h)
! solve y'+f(x,y) = 0 with 4th order Runge-Kutta
  real, intent(in) :: a, b, y0, h
  real :: x,y, k1, k2, k3, k4, e
  x=a; y=y0
  do while (x < b) 
    k1 = h * f(x,y)
    k2 = h * f(x+h/2, y+k1/2)
    k3 = h * f(x+h/2, y+k2/2)
    k4 = h * f(x+h, y+k3)
    y = y + (k1 + 2*k2 + 2*k3 + k4)/6.0
    x = x+h
  end do
  e=exp(-b)+b-1
  e=abs((y-e)/e)
  write(*,'(4F12.6)') h,x,y,e
end subroutine
end program

    0.100000    1.000000    0.367880    0.000001
    0.010000    1.009999    0.374219    0.017232
    0.001000    1.000991    0.368509    0.001712
    0.000100    1.000054    0.367884    0.000012

! modify the loop
  i=1; x=a; y=y0
  do while (x < b) 
    ...
    x = a+i*h
    i=i+1
  end do

    0.100000    1.000000    0.367880    0.000001
    0.010000    1.000000    0.367879    0.000000
    0.001000    1.000000    0.367880    0.000000
    0.000100    1.000000    0.367879    0.000000