module para
  implicit none
  integer :: i,j,k
  integer,parameter :: m = 21, &
                     & n = 10001
  real(8),parameter :: kappa = 1.d0, &
                     & tt = 10.d0,   &
                     & t0 = 0.d0,    &
                     & xx = 1.d0,    &
                     & x0 = 0.d0
  real(8) :: r,dx,dt
end module para
!*********************************************************************
program partial_equation
  use para
  implicit none
  real(8),dimension(0:m,0:n) :: u
  
  call calc_param
  call icon(u)
  call bcon(u)
  call calc_nextu(u)
  
  stop
  
end program partial_equation
!*******calculate parameter*******************************************
subroutine calc_param
  use para
  implicit none
  
  dt = (tt-t0)/dble(n-1)
  dx = (xx-x0)/dble(m-1)
  r  = kappa*dt/(dx*dx)
  
  return
end subroutine calc_param
!********initial condition********************************************
subroutine icon(u)
  use para
  implicit none
  real(8),dimension(0:m,0:n) :: u
  
  do i = 0,m
    if(i /= m) then
      u(i,0) = 0.d0
    else
      u(i,0) = 1.d0
    end if
  end do
  
  return
end subroutine icon
!*************boundary condition**************************************
subroutine bcon(u)
  use para
  implicit none
  real(8),dimension(0:m,0:n) :: u
  
  do j = 0,n
    u(0,j) = 0.d0
    u(m,j) = 1.d0
  end do
  
  return
end subroutine bcon
!*********************************************************************
subroutine calc_nextu(u)
  use para
  implicit none
  real(8),dimension(0:m,0:n) :: u
  
  do k = 1,n-1
    do j = 1,m-1
      u(j, k+1) = (1.0-2.0*r)*u(j,k)+r*(u(j+1,k)+u(j-1,k))
    end do
  end do
  
  return
  
end subroutine calc_nextu