Dim dx, p, cp, k, h, dt, Tin, Tout, i
Dim T(0 To 3), Tp(0 To 3)

'set properties
dx = 0.05 'm
p = 1860 'kg/m3
cp = 780 'J/kg K
k = 0.72 'W/m K
h = 10 'W/m2 K
dt = 1 's

'set init conditions
Tin = 20
Tout = 20
For i = 0 To 3
   T(i) = 20
Next i

'open output file
Open "c:\classes\mee410-00\finitediffwengystorage\temps.txt" For Output As #1

'solve for temps at end of time interval
For i = 1 To 200000
   'set new outside temp
   Tout = 0
   
   'find temp of each node at end of time interval
   Tp(0) = T(0) + (dt / (dx * p * cp)) * (h * (Tout - T(0)) + k * (T(1) - T(0)) / dx)
   Tp(1) = T(1) + (dt / (dx ^ 2 * p * cp)) * (k * (T(0) - T(1)) + k * (T(2) - T(1)))
   Tp(2) = T(2) + (dt / (dx ^ 2 * p * cp)) * (k * (T(1) - T(2)) + k * (T(3) - T(2)))
   Tp(3) = T(3) + (dt / (dx * p * cp)) * (h * (Tin - T(3)) + k * (T(2) - T(3)) / dx)
   
   'reset temps
   T(0) = Tp(0)
   T(1) = Tp(1)
   T(2) = Tp(2)
   T(3) = Tp(3)

   'print output to screen and file
   If i = 1000 Or (i = Int(i / 10000) * 10000) Then
      Picture1.Print i, Format(T(0), "0.00"), Format(T(1), "0.00"), Format(T(2), "0.00"), Format(T(3), "0.00")
      Print #1, i, Format(T(0), "0.00"), Format(T(1), "0.00"), Format(T(2), "0.00"), Format(T(3), "0.00")
   End If
Next i

Close