program heat

   integer :: j
   integer, parameter :: jbar=50
   double precision :: dx=dble(1.0), dt=dble(0.1), sig=dble(1.0), endtime=dble(10.0), time=dble(0.0)
   double precision :: TL = dble(400.0), TR = dble(0.0)
   double precision :: T(0:jbar+1), Tnew(0:jbar+1)

   T(:) = dble(0.0)

   do

      if ( time > endtime ) exit
      
      T(0)      = dble(2.0)*TL - T(1)     ! Boundary conditions
      T(jbar+1) = dble(2.0)*TR - T(jbar)

      do j=1,jbar
         Tnew(j) = T(j) + sig*dt*(T(j+1)+T(j-1)-dble(2.0)*T(j))/dx/dx
      enddo

      T(1:jbar) = Tnew(1:jbar)
      time = time + dt

   end do

   do j=1,jbar
      write(*,'(2(a,f11.6))') "x ",(dble(j-1)+dble(0.5))*dx," T ", T(j)
   end do

end program