1.
program adams
! solve y'=x-y0 using Adams-Bashforth method
! initial value y(0)=0
  call adamsb(0.0, 1.0,  0.05)

contains

real function exact(x)
  real, intent(in) :: x
  exact=exp(-x)+x-1.0
end function

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


subroutine adamsb(a, b, h)
  implicit none
  real, intent(in) :: a, b,  h
  real :: x, y, y1, y2, f0, f1, f2

  ! initial values
  f2 = f(a, exact(a))
  f1 = f(a+h, exact(a+h))
  x = a+2*h
  y = exact(x)
  f0 = f(x, y)

  do while (x <= b)
    ! predictor
    y1 = y + h/12*(23*f0-16*f1+5*f2) 

    x = x+h
    f2=f1; f1=f0; f0=f(x,y1)

    ! corrector
    y2 = y +h/12*(5*f0+8*f1-f2)
     
    y = y2

    write(*,'(3F10.6)') x,y,exact(x)
  end do
end subroutine

end program


  0.150000  0.010708  0.010708
  0.200000  0.018731  0.018731
  0.250000  0.028802  0.028801
  0.300000  0.040819  0.040818
  0.350000  0.054689  0.054688
  0.400000  0.070322  0.070320
  0.450000  0.087630  0.087628
  0.500000  0.106533  0.106531
  0.550000  0.126952  0.126950
  0.600000  0.148814  0.148812
  0.650000  0.172048  0.172046
  0.700000  0.196588  0.196585
  0.750000  0.222369  0.222367
  0.800000  0.249332  0.249329
  0.850000  0.277418  0.277415
  0.900000  0.306572  0.306570
  0.950000  0.336744  0.336741


2.
program shoot
! solve a boundary value problem using the shooting method
  implicit none
  integer, parameter :: nsize=10000
  real, parameter :: pi = 3.141592654
  real, dimension(0:nsize) :: x, y, v
  real :: y0, h, v0, v1, v2, yb, yb0, yb1
  integer :: i, nstep

  nstep=10
  h  = (pi/2)/nstep
  y0 = 1.0
  yb = 0.0

! y'(0) = 0
  v0=0.0
  call runge(y0, v0, nstep, h, x(0:nstep), y(0:nstep), v(0:nstep))
  write(*,'(/"v(0)=",F10.3)') v0
  do i=0,nstep
     write(*,'(3F10.3)') x(i), y(i),v(i)
  end do
  yb0=y(nstep)

! y'(0) = 1
  v1=1.0
  call runge(y0, v1, nstep, h, x(0:nstep), y(0:nstep), v(0:nstep))
  write(*,'(/"v(0)=",F10.3)') v1
  do i=0,nstep
     write(*,'(3F10.3)') x(i), y(i),v(i)
  end do
  yb1=y(nstep)

! y'(0) = v0
  v2 = v1-(yb1-yb)*(v1-v0)/(yb1-yb0)
  call runge(y0, v2, nstep, h, x(0:nstep), y(0:nstep), v(0:nstep))
  write(*,'(/"v(0)=",F10.3)') v2
  do i=0,nstep
     write(*,'(4F10.3)') x(i), y(i),v(i), exact(x(i))
  end do

contains

real function exact(x)
  real, intent(in):: x
  exact = x+cos(x)-pi/2*sin(x)
end function


real function f(x, y)
! u'+y=x
  real, intent(in) :: x, y
  f = x-y
end function

subroutine runge(y0, v0, n, h, xx, yy, vv)
! 4th order Runge-Kutta
! solve the pair of equations
!    y' = u,  u'=f(x, y)
! corresponding to y''= f(x, y) 
! initial values are y(0)=y0, y'(0)=v0
! the solution is returned in xx, yy, vv 
  real, intent(in) :: y0, v0
  real, intent(in) :: h
  integer, intent(in) :: n
  real, intent(out), dimension(0:n) :: xx, yy, vv
  integer i
  real :: x, y, v, ky1, ky2, ky3, ky4, kv1, kv2, kv3, kv4 

  x=0.0;  xx(0)=x 
  y=y0 ;  yy(0)=y0
  v=v0 ;  vv(0)=v0

  do i=1,n
    ky1 = h * v 
    kv1 = h * f(x, y) 

    ky2 = h * (v+kv1/2)
    kv2 = h * f(x+h/2, y+ky1/2)

    ky3 = h * (v+kv2/2)
    kv3 = h * f(x+h/2, y+ky2/2)

    ky4 = h * (v+kv3)
    kv4 = h * f(x+h, y+ky3)

    x = x+h
    y = y + (ky1 + 2*ky2 + 2*ky3 + ky4)/6.0
    v = v + (kv1 + 2*kv2 + 2*kv3 + kv4)/6.0 

    xx(i)=x
    yy(i)=y
    vv(i)=v
  end do

end subroutine

end program


v(0)=     0.000
     0.000     1.000     0.000
     0.157     0.988    -0.144
     0.314     0.956    -0.260
     0.471     0.908    -0.345
     0.628     0.850    -0.397
     0.785     0.785    -0.414
     0.942     0.721    -0.397
     1.100     0.663    -0.345
     1.257     0.615    -0.260
     1.414     0.582    -0.144
     1.571     0.571    -0.000

v(0)=     1.000
     0.000     1.000     1.000
     0.157     1.145     0.844
     0.314     1.265     0.691
     0.471     1.362     0.546
     0.628     1.437     0.412
     0.785     1.493     0.293
     0.942     1.530     0.191
     1.100     1.554     0.109
     1.257     1.566     0.049
     1.414     1.570     0.012
     1.571     1.571     0.000

v(0)=    -0.571
     0.000     1.000    -0.571     1.000
     0.157     0.899    -0.708     0.899
     0.314     0.780    -0.803     0.780
     0.471     0.649    -0.854     0.649
     0.628     0.514    -0.859     0.514
     0.785     0.382    -0.818     0.382
     0.942     0.259    -0.732     0.259
     1.100     0.154    -0.604     0.154
     1.257     0.072    -0.436     0.072
     1.414     0.019    -0.233     0.019
     1.571     0.000    -0.000     0.000


4.
program diff
! solve a boundary value problem using a difference method
  use linear
  implicit none

  integer, parameter :: n=11
  real, parameter :: pi=3.141592654
  real, dimension(n, n) :: A
  real, dimension(n) :: x, xx, y
  real :: h
  integer :: i

  h = (pi/2)/(n-1)

  do i=1,n
    x(i)=h*(i-1)
  end do

! coefficient matrix
  A=0.0
  A(1,1)=1.0
  do i=2,n-1
     A(i,i-1) = 1.0
     A(i,i)   = h**2-2.0
     A(i,i+1) = 1.0
  end do
  A(n,n)=1.0

! right hand side
  xx=x * h**2
  xx(1)=1.0
  xx(n)=0.0


! solve the set of linear equations

  call lu(A, n)
  call solve(A, xx, n)
  do i=1,n
    write(*,'(2F10.4)') x(i),xx(i)
  end do

end program


    0.0000    1.0000
    0.1571    0.8990
    0.3142    0.7797
    0.4712    0.6490
    0.6283    0.5138
    0.7854    0.3815
    0.9425    0.2591
    1.0996    0.1536
    1.2566    0.0714
    1.4137    0.0185
    1.5708    0.0000