! FTCS Method for solving Parabolic PDE
! Mustafa Khawari
! 17.10.99

implicit none
real, parameter :: dt = 0.1, a = 0.01
integer, parameter :: Imax = 21, Nmax = 500
real :: X (Imax), T (Imax, Nmax)
real :: tim, d, dx, XX
integer :: i, n, h

XX = 1.
dx = XX / (Imax-1)
h = int (Imax/2) + 1
d = (a * dt) / (dx**2)

do i = 1, Imax
  X (i) = (i-1) * dx
end do

do i = 1, Imax
  if (i <= h) then
    T (i, 1) = 200 * X (i)
  else
    T (i, 1) = 200 - 200 * X (i)
  end if
end do


open (6, file = "Para_1D_FTCS_EM.plt")
write (6, 10)

do n = 1, Nmax - 1
  tim = (n-1) * dt
  write (6, 11) tim
  
  do i = 1, Imax
    if (i == 1 .or. i == Imax) then
      T (i, n+1) = T (i, n)
    else
      T (i, n+1) = T (i, n) + d * (T (i+1, n) - 2 * T (i, n) + T (i-1, n))
    end if
    write (6, 12) X (i), T (i, n)
  end do
  
  print *, n
end do
print *, "     d =", d

10    format ( 'VARIABLES=X  C' )
11    format ( 'ZONE T="' , f5.1 , '"' )
12    format ( f10.4 , f100.4 )
end