Ecuación Del Calor En Fortran

Program ec_calor

interface
Subroutine DerSeg(U,F)
real*8,intent(in) :: U(:)
real*8,intent(inout) :: F(:)
endsubroutine
end interface

implicit none
real*8,parameter :: dt=0.01
integer,parameter :: n = 20,t = 100
real*8 ::U(1:2*n-2),F(1:2*n-2),Fm1(1:2*n-2),Fp(1:2*n-2),x(0:n),g,U0,Un
integer :: i,jopen(30,file="PRCT9VAL.dat",status="unknown",action="write")
!----CC,CI
U0 = 0
Un = 0
U = 0
F = 0
Fm1= 0
x(0)=0
x(n)=1
do i = 1,n-1
x(i) = 1d0*i/nU(i) = g(x(i))
end do

call DerSeg(U,Fm1)

!----Euler primer paso
U = U + dt*Fm1
F(1:n-1) = U(n:2*(n-1))
callDerSeg(U,F)

!----AdamsBaschford
write(30,"(110F17.3)")0d0,x(:)
do i = 1,t
write(*,*) "Estado Bucle",i
!U predecidaU = U + (3d0*F-Fm1)*dt/2d0
!Calculo F predecida
Fp(1:n-1) = U(n:2*(n-1))
call DerSeg(U,Fp)
!Calculo U Siguiente
U =U + (F + Fp)*dt/2d0
!Calculo nueva y antigua F
Fm1 = F
F(1:n-1) = U(n:2*(n-1))
call DerSeg(U,F)!Escribo vector para t en el que estoy

write(30,"(110F17.6)")dt*i,U0,U(1:n-1),Un

end do
end program



Function g(x)real*8,intent(in) :: x
real*8 :: g

!g=exp((-((x-1d0/2)/0.1d0)**2))

g=0
if ( x >= 0.25 .and. x 0.5 .and. x