program transient_fdm_solver
    implicit none

    ! Variable declarations
    integer :: geom_type, n_nodes, method
    integer :: bc_left_type, bc_right_type, init_type
    real(8) :: length, dt_choice, t_sim
    real(8) :: k_0, rho, cp_0, beta_k, beta_cp
    real(8) :: bc_left_p1, bc_left_p2, bc_right_p1, bc_right_p2
    real(8) :: init_p1, init_p2

    ! Grid and solver variables
    real(8), allocatable, dimension(:) :: T, T_old, r_nodes, V_nodes
    real(8), allocatable, dimension(:) :: a_coeff, b_coeff, c_coeff, d_coeff
    real(8), allocatable, dimension(:) :: k_nodes, cp_nodes
    real(8) :: dx, dt, alpha, current_time
    integer :: i, step, total_steps, print_interval
    character(len=20) :: geom_name, method_name
    real(8) :: E_in_total, E_out_total, E_st_initial, E_st_current, E_st_change
    real(8) :: balance_err, heat_flux_left, heat_flux_right
    real(8) :: dt_max_stable, Fo_max, Fo_node, Bi_node
    logical :: is_stable

    ! Read parameters from stdin
    read *, geom_type      ! 1=Plane, 2=Cylinder, 3=Sphere
    read *, length         ! Total thickness L (m) or outer radius r0 (m)
    read *, n_nodes        ! Number of nodes
    read *, dt_choice      ! <= 0 for auto stability, > 0 for manual time step (s)
    read *, method         ! 1=Explicit, 2=Implicit BE, 3=Crank-Nicolson
    read *, k_0            ! Thermal conductivity k0 [W/m-K]
    read *, rho            ! Density rho [kg/m3]
    read *, cp_0           ! Specific heat cp0 [J/kg-K]
    read *, beta_k         ! Temperature coeff for conductivity
    read *, beta_cp        ! Temperature coeff for specific heat
    read *, bc_left_type   ! Left BC: 1=Const T, 2=Const q", 3=Convection, 4=Insulated
    read *, bc_left_p1     ! T_s, q", or h
    read *, bc_left_p2     ! T_inf (if convection)
    read *, bc_right_type  ! Right BC: 1=Const T, 2=Const q", 3=Convection, 4=Insulated
    read *, bc_right_p1    ! T_s, q", or h
    read *, bc_right_p2    ! T_inf (if convection)
    read *, init_type      ! 1=Uniform, 2=Linear, 3=Custom (parabolic)
    read *, init_p1        ! Ti or T_left or T_center
    read *, init_p2        ! T_right or T_outer
    read *, t_sim          ! Total simulation time [s]

    ! Allocations
    allocate(T(n_nodes))
    allocate(T_old(n_nodes))
    allocate(r_nodes(n_nodes))
    allocate(V_nodes(n_nodes))
    allocate(a_coeff(n_nodes))
    allocate(b_coeff(n_nodes))
    allocate(c_coeff(n_nodes))
    allocate(d_coeff(n_nodes))
    allocate(k_nodes(n_nodes))
    allocate(cp_nodes(n_nodes))

    ! Grid setup
    dx = length / dble(n_nodes - 1)
    do i = 1, n_nodes
        r_nodes(i) = dble(i - 1) * dx
    end do

    ! Calculate volumes for nodes (for energy balance check)
    select case (geom_type)
        case (1)
            geom_name = "1D Plane Wall"
            do i = 1, n_nodes
                if (i == 1 .or. i == n_nodes) then
                    V_nodes(i) = dx / 2.0d0
                else
                    V_nodes(i) = dx
                end if
            end do
        case (2)
            geom_name = "1D Cylinder"
            do i = 1, n_nodes
                if (i == 1) then
                    V_nodes(i) = 3.14159265358979d0 * (dx / 2.0d0)**2
                elseif (i == n_nodes) then
                    V_nodes(i) = 3.14159265358979d0 * (length**2 - (length - dx/2.0d0)**2)
                else
                    V_nodes(i) = 3.14159265358979d0 * ((r_nodes(i) + dx/2.0d0)**2 - (r_nodes(i) - dx/2.0d0)**2)
                end if
            end do
        case (3)
            geom_name = "1D Sphere"
            do i = 1, n_nodes
                if (i == 1) then
                    V_nodes(i) = (4.0d0/3.0d0) * 3.14159265358979d0 * (dx / 2.0d0)**3
                elseif (i == n_nodes) then
                    V_nodes(i) = (4.0d0/3.0d0) * 3.14159265358979d0 * (length**3 - (length - dx/2.0d0)**3)
                else
                    V_nodes(i) = (4.0d0/3.0d0) * 3.14159265358979d0 * ((r_nodes(i) + dx/2.0d0)**3 - (r_nodes(i) - dx/2.0d0)**3)
                end if
            end do
    end select

    ! Initial temperature distribution
    do i = 1, n_nodes
        select case (init_type)
            case (1) ! Uniform
                T(i) = init_p1
            case (2) ! Linear
                T(i) = init_p1 + (init_p2 - init_p1) * (r_nodes(i) / length)
            case (3) ! Parabolic
                T(i) = init_p1 - (init_p1 - init_p2) * (r_nodes(i) / length)**2
        end select
    end do

    ! Estimate stability limit for Explicit method
    ! We evaluate at initial temperature (which is safe)
    call get_properties(T(1), k_0, cp_0, beta_k, beta_cp, k_nodes(1), cp_nodes(1))
    alpha = k_nodes(1) / (rho * cp_nodes(1))
    select case (geom_type)
        case (1)
            dt_max_stable = (dx**2) / (2.0d0 * alpha)
            ! Adjust for convection if present
            if (bc_left_type == 3) then
                dt_max_stable = min(dt_max_stable, (dx**2) / (2.0d0 * alpha * (1.0d0 + bc_left_p1 * dx / k_nodes(1))))
            endif
            if (bc_right_type == 3) then
                dt_max_stable = min(dt_max_stable, (dx**2) / (2.0d0 * alpha * (1.0d0 + bc_right_p1 * dx / k_nodes(1))))
            endif
        case (2)
            ! Cylinder center node has 4*Fo limit
            dt_max_stable = (dx**2) / (4.0d0 * alpha)
        case (3)
            ! Sphere center node has 6*Fo limit
            dt_max_stable = (dx**2) / (6.0d0 * alpha)
    end select

    ! Time step selection
    if (dt_choice <= 0.0d0) then
        dt = 0.95d0 * dt_max_stable
    else
        dt = dt_choice
    end if

    ! Determine method name
    select case (method)
        case (1)
            method_name = "Explicit (Euler)"
        case (2)
            method_name = "Implicit (Backward)"
        case (3)
            method_name = "Crank-Nicolson"
    end select

    total_steps = ceiling(t_sim / dt)
    if (total_steps < 1) total_steps = 1
    
    ! Adjust print interval to output ~10 profiles in text report
    print_interval = max(1, total_steps / 10)

    ! Initial energy stored
    E_st_initial = 0.0d0
    do i = 1, n_nodes
        call get_properties(T(i), k_0, cp_0, beta_k, beta_cp, k_nodes(i), cp_nodes(i))
        E_st_initial = E_st_initial + rho * cp_nodes(i) * T(i) * V_nodes(i)
    end do

    E_in_total = 0.0d0
    E_out_total = 0.0d0
    current_time = 0.0d0

    ! ==========================================
    ! OUTPUT HEADER REPORT
    ! ==========================================
    print *, '==================================================='
    print *, '  1D TRANSIENT CONDUCTION NUMERICAL SOLVER REPORT'
    print *, '==================================================='
    print *, ''
    print '(A,A)', '  Geometry Configuration:      ', trim(geom_name)
    print '(A,F10.6,A)', '  Domain Length/Radius (L):   ', length, ' m'
    print '(A,I5)', '  Number of Nodal Points (N):  ', n_nodes
    print '(A,F10.6,A)', '  Nodal Grid Spacing (dx):    ', dx, ' m'
    print '(A,A)', '  Time-Marching Algorithm:     ', trim(method_name)
    print '(A,F12.6,A)', '  Time Step Size (dt):        ', dt, ' s'
    print '(A,I8)', '  Total Number of Steps:       ', total_steps
    print '(A,F12.4,A)', '  Total Simulation Time:      ', t_sim, ' s'
    print *, ''
    
    print *, '==================================================='
    print *, '  MATERIAL PROPERTIES & BOUNDARY CONDITIONS'
    print *, '==================================================='
    print '(A,F10.4,A)', '  Thermal Conductivity k0:    ', k_0, ' W/(m-K)'
    if (beta_k /= 0.0d0) then
        print '(A,ES12.4,A)', '    Temperature Coeff (beta_k):', beta_k, ' 1/K'
    endif
    print '(A,F10.2,A)', '  Density (rho):              ', rho, ' kg/m3'
    print '(A,F10.2,A)', '  Specific Heat Cp0:          ', cp_0, ' J/(kg-K)'
    if (beta_cp /= 0.0d0) then
        print '(A,ES12.4,A)', '    Temperature Coeff (beta_cp):', beta_cp, ' 1/K'
    endif
    print *, ''
    
    ! Boundary Condition displays
    call print_bc_desc("Left/Inner (x=0)", bc_left_type, bc_left_p1, bc_left_p2)
    call print_bc_desc("Right/Outer (x=L)", bc_right_type, bc_right_p1, bc_right_p2)
    print *, ''

    ! Stability indicator
    is_stable = .true.
    if (method == 1 .and. dt > dt_max_stable) then
        is_stable = .false.
        print *, '⚠️ WARNING: Explicit stability limit exceeded!'
        print '(A,F12.6,A)', '  Calculated stable dt limit is ', dt_max_stable, ' s'
        print '(A,F12.6,A)', '  Your selected time step dt is ', dt, ' s'
    else
        print '(A,F12.6,A)', '  Maximum stable explicit dt limit: ', dt_max_stable, ' s'
        print *, '  Solver stability status: OK'
    endif
    print *, ''

    print *, '==================================================='
    print *, '  TRANSIENT TEMPERATURE FIELD PROGRESSION'
    print *, '==================================================='
    print '(A,9A10)', ' Time [s]', 'Node 1', 'Node 25%', 'Node 50%', 'Node 75%', 'Node N', 'E_in[J]', 'E_st[J]', 'Err [%]'
    
    ! Print initial state
    call compute_energy_change(T, T_old, V_nodes, rho, k_0, cp_0, beta_k, beta_cp, &
                               bc_left_type, bc_left_p1, bc_left_p2, &
                               bc_right_type, bc_right_p1, bc_right_p2, &
                               0.0d0, 0.0d0, E_st_initial, E_st_current, &
                               E_in_total, E_out_total, balance_err)
    call print_state_line(0.0d0, T, n_nodes, E_in_total, E_out_total, E_st_current, E_st_initial, balance_err)

    ! Time Loop
    do step = 1, total_steps
        T_old = T
        current_time = dble(step) * dt

        ! Update node properties
        do i = 1, n_nodes
            call get_properties(T_old(i), k_0, cp_0, beta_k, beta_cp, k_nodes(i), cp_nodes(i))
        end do

        select case (method)
            case (1) ! Explicit Euler
                call solve_explicit(T, T_old, n_nodes, geom_type, dx, dt, rho, k_nodes, cp_nodes, &
                                    bc_left_type, bc_left_p1, bc_left_p2, &
                                    bc_right_type, bc_right_p1, bc_right_p2)
            case (2, 3) ! Implicit BE or Crank-Nicolson
                call solve_implicit(T, T_old, n_nodes, geom_type, dx, dt, rho, k_nodes, cp_nodes, method, &
                                    bc_left_type, bc_left_p1, bc_left_p2, &
                                    bc_right_type, bc_right_p1, bc_right_p2)
        end select

        ! Compute border heat fluxes and cumulative energy inputs (for 1 m2 surface)
        call get_boundary_fluxes(T, T_old, n_nodes, geom_type, dx, k_nodes(1), k_nodes(n_nodes), &
                                 bc_left_type, bc_left_p1, bc_left_p2, &
                                 bc_right_type, bc_right_p1, bc_right_p2, &
                                 heat_flux_left, heat_flux_right)
        
        E_in_total = E_in_total + heat_flux_left * dt
        E_out_total = E_out_total + heat_flux_right * dt

        if (mod(step, print_interval) == 0 .or. step == total_steps) then
            call compute_energy_change(T, T_old, V_nodes, rho, k_0, cp_0, beta_k, beta_cp, &
                                       bc_left_type, bc_left_p1, bc_left_p2, &
                                       bc_right_type, bc_right_p1, bc_right_p2, &
                                       E_in_total, E_out_total, E_st_initial, E_st_current, &
                                       E_in_total, E_out_total, balance_err)
            call print_state_line(current_time, T, n_nodes, E_in_total, E_out_total, E_st_current, E_st_initial, balance_err)
        end if
    end do

    print *, ''
    print *, '==================================================='
    print *, '  FINAL NODAL TEMPERATURE DISTRIBUTION'
    print *, '==================================================='
    print '(A12,A18,A20)', ' Node Index', 'Coordinate [m]', 'Temperature [C]'
    do i = 1, n_nodes
        if (i == 1 .or. i == n_nodes .or. mod(i, max(1, n_nodes/20)) == 0) then
            print '(I10,F18.6,F20.4)', i, r_nodes(i), T(i)
        end if
    end do
    print *, '==================================================='

contains

    subroutine get_properties(Temp, k0, cp0, beta_k, beta_cp, k_out, cp_out)
        real(8), intent(in) :: Temp, k0, cp0, beta_k, beta_cp
        real(8), intent(out) :: k_out, cp_out
        k_out = k0 * (1.0d0 + beta_k * (Temp - 20.0d0))
        cp_out = cp0 * (1.0d0 + beta_cp * (Temp - 20.0d0))
    end subroutine get_properties

    subroutine print_bc_desc(side, bc_type, p1, p2)
        character(len=*), intent(in) :: side
        integer, intent(in) :: bc_type
        real(8), intent(in) :: p1, p2
        select case (bc_type)
            case (1)
                print '(A,A,F10.2,A)', '  BC at ', side, p1, ' deg-C (Constant Temperature)'
            case (2)
                print '(A,A,F10.2,A)', '  BC at ', side, p1, ' W/m2 (Constant Heat Flux)'
            case (3)
                print '(A,A,F10.2,A,F10.2,A)', '  BC at ', side, p1, ' W/m2-K (Convection), T_inf = ', p2, ' deg-C'
            case (4)
                print '(A,A,A)', '  BC at ', side, ' Insulated (Adiabatic)'
        end select
    end subroutine print_bc_desc

    subroutine get_boundary_fluxes(T_now, T_prev, n, geom, dx_val, k_left, k_right, &
                                   bc_l_type, bc_l_p1, bc_l_p2, &
                                   bc_r_type, bc_r_p1, bc_r_p2, &
                                   q_left, q_right)
        real(8), dimension(n), intent(in) :: T_now, T_prev
        integer, intent(in) :: n, geom, bc_l_type, bc_r_type
        real(8), intent(in) :: dx_val, k_left, k_right, bc_l_p1, bc_l_p2, bc_r_p1, bc_r_p2
        real(8), intent(out) :: q_left, q_right

        ! Left Boundary Inward Heat Flux
        select case (bc_l_type)
            case (1) ! Constant T
                q_left = k_left * (T_now(1) - T_now(2)) / dx_val
            case (2) ! Constant q"
                q_left = bc_l_p1
            case (3) ! Convection
                q_left = bc_l_p1 * (bc_l_p2 - T_now(1))
            case (4) ! Insulated
                q_left = 0.0d0
        end select

        ! Center symmetry check for radial domains
        if (geom > 1) q_left = 0.0d0

        ! Right Boundary Outward Heat Flux (positive escaping)
        select case (bc_r_type)
            case (1) ! Constant T
                q_right = k_right * (T_now(n) - T_now(n-1)) / dx_val
            case (2) ! Constant q" (inward flux so escaping is -inward)
                q_right = -bc_r_p1
            case (3) ! Convection
                q_right = bc_r_p1 * (T_now(n) - bc_r_p2)
            case (4) ! Insulated
                q_right = 0.0d0
        end select
    end subroutine get_boundary_fluxes

    subroutine print_state_line(t_val, T_arr, n, E_in, E_out, E_st, E_st_init, err)
        real(8), intent(in) :: t_val, E_in, E_out, E_st, E_st_init, err
        integer, intent(in) :: n
        real(8), dimension(n), intent(in) :: T_arr
        print '(F10.3,5F10.3,3ES10.2)', t_val, T_arr(1), T_arr(max(1,n/4)), T_arr(max(1,n/2)), &
                                         T_arr(max(1,3*n/4)), T_arr(n), (E_in - E_out), (E_st - E_st_init), err
    end subroutine print_state_line

    subroutine compute_energy_change(T_now, T_prev, V_arr, rho_val, k0, cp0, beta_k, beta_cp, &
                                     bc_l_type, bc_l_p1, bc_l_p2, &
                                     bc_r_type, bc_r_p1, bc_r_p2, &
                                     E_in, E_out, E_st_init, E_st, &
                                     E_in_out, E_out_out, err)
        real(8), dimension(:), intent(in) :: T_now, T_prev, V_arr
        real(8), intent(in) :: rho_val, k0, cp0, beta_k, beta_cp
        integer, intent(in) :: bc_l_type, bc_r_type
        real(8), intent(in) :: bc_l_p1, bc_l_p2, bc_r_p1, bc_r_p2, E_in, E_out, E_st_init
        real(8), intent(out) :: E_st, E_in_out, E_out_out, err

        integer :: idx
        real(8) :: temp_k, temp_cp, dE_net

        E_st = 0.0d0
        do idx = 1, size(T_now)
            call get_properties(T_now(idx), k0, cp0, beta_k, beta_cp, temp_k, temp_cp)
            E_st = E_st + rho_val * temp_cp * T_now(idx) * V_arr(idx)
        end do

        dE_net = E_in - E_out
        E_in_out = E_in
        E_out_out = E_out
        
        err = abs(dE_net - (E_st - E_st_init)) / (abs(E_st - E_st_init) + 1.0d0) * 100.0d0
    end subroutine compute_energy_change

    subroutine solve_explicit(T_new, T_now, n, geom, dx_val, dt_val, rho_val, k_arr, cp_arr, &
                              bc_l_type, bc_l_p1, bc_l_p2, &
                              bc_r_type, bc_r_p1, bc_r_p2)
        integer, intent(in) :: n, geom, bc_l_type, bc_r_type
        real(8), dimension(n), intent(in) :: T_now, k_arr, cp_arr
        real(8), dimension(n), intent(out) :: T_new
        real(8), intent(in) :: dx_val, dt_val, rho_val, bc_l_p1, bc_l_p2, bc_r_p1, bc_r_p2

        integer :: idx
        real(8) :: r_i, r_left, r_right, k_left, k_right, Fo_local

        ! Interior nodes update
        do idx = 2, n-1
            k_left = 0.5d0 * (k_arr(idx-1) + k_arr(idx))
            k_right = 0.5d0 * (k_arr(idx) + k_arr(idx+1))
            
            select case (geom)
                case (1) ! Plane
                    T_new(idx) = T_now(idx) + (dt_val / (rho_val * cp_arr(idx) * dx_val**2)) * &
                                 (k_left * (T_now(idx-1) - T_now(idx)) + k_right * (T_now(idx+1) - T_now(idx)))
                case (2) ! Cylinder
                    r_i = dble(idx - 1) * dx_val
                    r_left = r_i - 0.5d0 * dx_val
                    r_right = r_i + 0.5d0 * dx_val
                    T_new(idx) = T_now(idx) + (dt_val / (rho_val * cp_arr(idx) * r_i * dx_val**2)) * &
                                 (r_left * k_left * (T_now(idx-1) - T_now(idx)) + r_right * k_right * (T_now(idx+1) - T_now(idx)))
                case (3) ! Sphere
                    r_i = dble(idx - 1) * dx_val
                    r_left = r_i - 0.5d0 * dx_val
                    r_right = r_i + 0.5d0 * dx_val
                    T_new(idx) = T_now(idx) + (dt_val / (rho_val * cp_arr(idx) * r_i**2 * dx_val**2)) * &
                                 (r_left**2 * k_left * (T_now(idx-1) - T_now(idx)) + r_right**2 * k_right * (T_now(idx+1) - T_now(idx)))
            end select
        end do

        ! Left/Inner boundary node (idx=1)
        if (geom > 1) then
            ! Center node (insulated symmetry boundary)
            Fo_local = k_arr(1) * dt_val / (rho_val * cp_arr(1) * dx_val**2)
            if (geom == 2) then
                T_new(1) = T_now(1) + 4.0d0 * Fo_local * (T_now(2) - T_now(1))
            else ! Sphere
                T_new(1) = T_now(1) + 6.0d0 * Fo_local * (T_now(2) - T_now(1))
            end if
        else
            select case (bc_l_type)
                case (1) ! Constant T
                    T_new(1) = bc_l_p1
                case (2) ! Constant q"
                    Fo_local = k_arr(1) * dt_val / (rho_val * cp_arr(1) * dx_val**2)
                    T_new(1) = T_now(1) + 2.0d0 * Fo_local * (T_now(2) - T_now(1) + bc_l_p1 * dx_val / k_arr(1))
                case (3) ! Convection
                    Fo_local = k_arr(1) * dt_val / (rho_val * cp_arr(1) * dx_val**2)
                    T_new(1) = T_now(1) + 2.0d0 * Fo_local * (T_now(2) - T_now(1) + bc_l_p1 * dx_val * (bc_l_p2 - T_now(1)) / k_arr(1))
                case (4) ! Insulated
                    Fo_local = k_arr(1) * dt_val / (rho_val * cp_arr(1) * dx_val**2)
                    T_new(1) = T_now(1) + 2.0d0 * Fo_local * (T_now(2) - T_now(1))
            end select
        end if

        ! Right/Outer boundary node (idx=n)
        select case (bc_r_type)
            case (1) ! Constant T
                T_new(n) = bc_r_p1
            case (2) ! Constant q"
                Fo_local = k_arr(n) * dt_val / (rho_val * cp_arr(n) * dx_val**2)
                select case (geom)
                    case (1)
                        T_new(n) = T_now(n) + 2.0d0 * Fo_local * (T_now(n-1) - T_now(n) + bc_r_p1 * dx_val / k_arr(n))
                    case (2)
                        T_new(n) = T_now(n) + 2.0d0 * Fo_local * ((1.0d0 - dx_val/(2.0d0*length)) * &
                                   (T_now(n-1) - T_now(n)) + bc_r_p1 * dx_val / k_arr(n))
                    case (3)
                        T_new(n) = T_now(n) + 2.0d0 * Fo_local * ((1.0d0 - dx_val/length)**2 * &
                                   (T_now(n-1) - T_now(n)) + bc_r_p1 * dx_val / k_arr(n))
                end select
            case (3) ! Convection
                Fo_local = k_arr(n) * dt_val / (rho_val * cp_arr(n) * dx_val**2)
                select case (geom)
                    case (1)
                        T_new(n) = T_now(n) + 2.0d0 * Fo_local * (T_now(n-1) - T_now(n) + bc_r_p1 * dx_val * (bc_r_p2 - T_now(n)) / k_arr(n))
                    case (2)
                        T_new(n) = T_now(n) + 2.0d0 * Fo_local * ((1.0d0 - dx_val/(2.0d0*length)) * &
                                   (T_now(n-1) - T_now(n)) + bc_r_p1 * dx_val * (bc_r_p2 - T_now(n)) / k_arr(n))
                    case (3)
                        T_new(n) = T_now(n) + 2.0d0 * Fo_local * ((1.0d0 - dx_val/length)**2 * &
                                   (T_now(n-1) - T_now(n)) + bc_r_p1 * dx_val * (bc_r_p2 - T_now(n)) / k_arr(n))
                end select
            case (4) ! Insulated
                Fo_local = k_arr(n) * dt_val / (rho_val * cp_arr(n) * dx_val**2)
                select case (geom)
                    case (1)
                        T_new(n) = T_now(n) + 2.0d0 * Fo_local * (T_now(n-1) - T_now(n))
                    case (2)
                        T_new(n) = T_now(n) + 2.0d0 * Fo_local * (1.0d0 - dx_val/(2.0d0*length)) * (T_now(n-1) - T_now(n))
                    case (3)
                        T_new(n) = T_now(n) + 2.0d0 * Fo_local * (1.0d0 - dx_val/length)**2 * (T_now(n-1) - T_now(n))
                end select
        end select
    end subroutine solve_explicit

    subroutine solve_implicit(T_new, T_now, n, geom, dx_val, dt_val, rho_val, k_arr, cp_arr, meth_idx, &
                              bc_l_type, bc_l_p1, bc_l_p2, &
                              bc_r_type, bc_r_p1, bc_r_p2)
        integer, intent(in) :: n, geom, meth_idx, bc_l_type, bc_r_type
        real(8), dimension(n), intent(in) :: T_now, k_arr, cp_arr
        real(8), dimension(n), intent(out) :: T_new
        real(8), intent(in) :: dx_val, dt_val, rho_val, bc_l_p1, bc_l_p2, bc_r_p1, bc_r_p2

        real(8), dimension(n) :: a, b, c, d
        integer :: idx
        real(8) :: r_i, r_left, r_right, k_left, k_right, Fo_local
        real(8) :: factor, Bi_local

        factor = 1.0d0
        if (meth_idx == 3) factor = 0.5d0 ! Crank-Nicolson uses 0.5 weight for both implicit/explicit

        ! Interior nodes coefficients
        do idx = 2, n-1
            k_left = 0.5d0 * (k_arr(idx-1) + k_arr(idx))
            k_right = 0.5d0 * (k_arr(idx) + k_arr(idx+1))
            
            select case (geom)
                case (1) ! Plane
                    a(idx) = -factor * (dt_val * k_left / (rho_val * cp_arr(idx) * dx_val**2))
                    c(idx) = -factor * (dt_val * k_right / (rho_val * cp_arr(idx) * dx_val**2))
                case (2) ! Cylinder
                    r_i = dble(idx - 1) * dx_val
                    r_left = r_i - 0.5d0 * dx_val
                    r_right = r_i + 0.5d0 * dx_val
                    a(idx) = -factor * (dt_val * r_left * k_left / (rho_val * cp_arr(idx) * r_i * dx_val**2))
                    c(idx) = -factor * (dt_val * r_right * k_right / (rho_val * cp_arr(idx) * r_i * dx_val**2))
                case (3) ! Sphere
                    r_i = dble(idx - 1) * dx_val
                    r_left = r_i - 0.5d0 * dx_val
                    r_right = r_i + 0.5d0 * dx_val
                    a(idx) = -factor * (dt_val * r_left**2 * k_left / (rho_val * cp_arr(idx) * r_i**2 * dx_val**2))
                    c(idx) = -factor * (dt_val * r_right**2 * k_right / (rho_val * cp_arr(idx) * r_i**2 * dx_val**2))
            end select
            b(idx) = 1.0d0 - a(idx) - c(idx)
            
            if (meth_idx == 2) then
                d(idx) = T_now(idx)
            else ! Crank-Nicolson RHS includes explicit part
                d(idx) = T_now(idx) - a(idx) * T_now(idx-1) + (a(idx) + c(idx)) * T_now(idx) - c(idx) * T_now(idx+1)
            end if
        end do

        ! Left boundary (idx=1)
        a(1) = 0.0d0
        if (geom > 1) then
            ! Symmetry at center
            Fo_local = k_arr(1) * dt_val / (rho_val * cp_arr(1) * dx_val**2)
            if (geom == 2) then
                b(1) = 1.0d0 + 4.0d0 * factor * Fo_local
                c(1) = -4.0d0 * factor * Fo_local
            else ! Sphere
                b(1) = 1.0d0 + 6.0d0 * factor * Fo_local
                c(1) = -6.0d0 * factor * Fo_local
            end if
            if (meth_idx == 2) then
                d(1) = T_now(1)
            else
                d(1) = T_now(1) - c(1) * (T_now(2) - T_now(1))
            end if
        else
            select case (bc_l_type)
                case (1) ! Constant T
                    b(1) = 1.0d0
                    c(1) = 0.0d0
                    d(1) = bc_l_p1
                case (2) ! Constant q"
                    Fo_local = k_arr(1) * dt_val / (rho_val * cp_arr(1) * dx_val**2)
                    b(1) = 1.0d0 + 2.0d0 * factor * Fo_local
                    c(1) = -2.0d0 * factor * Fo_local
                    if (meth_idx == 2) then
                        d(1) = T_now(1) + 2.0d0 * Fo_local * bc_l_p1 * dx_val / k_arr(1)
                    else
                        d(1) = T_now(1) - c(1) * (T_now(2) - T_now(1)) + 2.0d0 * Fo_local * bc_l_p1 * dx_val / k_arr(1)
                    end if
                case (3) ! Convection
                    Fo_local = k_arr(1) * dt_val / (rho_val * cp_arr(1) * dx_val**2)
                    Bi_local = bc_l_p1 * dx_val / k_arr(1)
                    b(1) = 1.0d0 + 2.0d0 * factor * Fo_local * (1.0d0 + Bi_local)
                    c(1) = -2.0d0 * factor * Fo_local
                    if (meth_idx == 2) then
                        d(1) = T_now(1) + 2.0d0 * Fo_local * Bi_local * bc_l_p2
                    else
                        d(1) = T_now(1) - c(1) * (T_now(2) - T_now(1)) + &
                               Fo_local * Bi_local * (2.0d0 * bc_l_p2 - T_now(1))
                    end if
                case (4) ! Insulated
                    Fo_local = k_arr(1) * dt_val / (rho_val * cp_arr(1) * dx_val**2)
                    b(1) = 1.0d0 + 2.0d0 * factor * Fo_local
                    c(1) = -2.0d0 * factor * Fo_local
                    if (meth_idx == 2) then
                        d(1) = T_now(1)
                    else
                        d(1) = T_now(1) - c(1) * (T_now(2) - T_now(1))
                    end if
            end select
        end if

        ! Right boundary (idx=n)
        c(n) = 0.0d0
        select case (bc_r_type)
            case (1) ! Constant T
                a(n) = 0.0d0
                b(n) = 1.0d0
                d(n) = bc_r_p1
            case (2) ! Constant q"
                Fo_local = k_arr(n) * dt_val / (rho_val * cp_arr(n) * dx_val**2)
                select case (geom)
                    case (1)
                        a(n) = -2.0d0 * factor * Fo_local
                        b(n) = 1.0d0 + 2.0d0 * factor * Fo_local
                        if (meth_idx == 2) then
                            d(n) = T_now(n) + 2.0d0 * Fo_local * bc_r_p1 * dx_val / k_arr(n)
                        else
                            d(n) = T_now(n) - a(n) * (T_now(n-1) - T_now(n)) + 2.0d0 * Fo_local * bc_r_p1 * dx_val / k_arr(n)
                        end if
                    case (2)
                        a(n) = -2.0d0 * factor * Fo_local * (1.0d0 - dx_val/(2.0d0*length))
                        b(n) = 1.0d0 - a(n)
                        if (meth_idx == 2) then
                            d(n) = T_now(n) + 2.0d0 * Fo_local * bc_r_p1 * dx_val / k_arr(n)
                        else
                            d(n) = T_now(n) - a(n) * (T_now(n-1) - T_now(n)) + 2.0d0 * Fo_local * bc_r_p1 * dx_val / k_arr(n)
                        end if
                    case (3)
                        a(n) = -2.0d0 * factor * Fo_local * (1.0d0 - dx_val/length)**2
                        b(n) = 1.0d0 - a(n)
                        if (meth_idx == 2) then
                            d(n) = T_now(n) + 2.0d0 * Fo_local * bc_r_p1 * dx_val / k_arr(n)
                        else
                            d(n) = T_now(n) - a(n) * (T_now(n-1) - T_now(n)) + 2.0d0 * Fo_local * bc_r_p1 * dx_val / k_arr(n)
                        end if
                end select
            case (3) ! Convection
                Fo_local = k_arr(n) * dt_val / (rho_val * cp_arr(n) * dx_val**2)
                Bi_local = bc_r_p1 * dx_val / k_arr(n)
                select case (geom)
                    case (1)
                        a(n) = -2.0d0 * factor * Fo_local
                        b(n) = 1.0d0 + 2.0d0 * factor * Fo_local * (1.0d0 + Bi_local)
                        if (meth_idx == 2) then
                            d(n) = T_now(n) + 2.0d0 * Fo_local * Bi_local * bc_r_p2
                        else
                            d(n) = T_now(n) - a(n) * (T_now(n-1) - T_now(n)) + &
                                   Fo_local * Bi_local * (2.0d0 * bc_r_p2 - T_now(n))
                        end if
                    case (2)
                        a(n) = -2.0d0 * factor * Fo_local * (1.0d0 - dx_val/(2.0d0*length))
                        b(n) = 1.0d0 - a(n) + 2.0d0 * factor * Fo_local * Bi_local
                        if (meth_idx == 2) then
                            d(n) = T_now(n) + 2.0d0 * Fo_local * Bi_local * bc_r_p2
                        else
                            d(n) = T_now(n) - a(n) * (T_now(n-1) - T_now(n)) + &
                                   Fo_local * Bi_local * (2.0d0 * bc_r_p2 - T_now(n))
                        end if
                    case (3)
                        a(n) = -2.0d0 * factor * Fo_local * (1.0d0 - dx_val/length)**2
                        b(n) = 1.0d0 - a(n) + 2.0d0 * factor * Fo_local * Bi_local
                        if (meth_idx == 2) then
                            d(n) = T_now(n) + 2.0d0 * Fo_local * Bi_local * bc_r_p2
                        else
                            d(n) = T_now(n) - a(n) * (T_now(n-1) - T_now(n)) + &
                                   Fo_local * Bi_local * (2.0d0 * bc_r_p2 - T_now(n))
                        end if
                end select
            case (4) ! Insulated
                Fo_local = k_arr(n) * dt_val / (rho_val * cp_arr(n) * dx_val**2)
                select case (geom)
                    case (1)
                        a(n) = -2.0d0 * factor * Fo_local
                        b(n) = 1.0d0 + 2.0d0 * factor * Fo_local
                    case (2)
                        a(n) = -2.0d0 * factor * Fo_local * (1.0d0 - dx_val/(2.0d0*length))
                        b(n) = 1.0d0 - a(n)
                    case (3)
                        a(n) = -2.0d0 * factor * Fo_local * (1.0d0 - dx_val/length)**2
                        b(n) = 1.0d0 - a(n)
                end select
                if (meth_idx == 2) then
                    d(n) = T_now(n)
                else
                    d(n) = T_now(n) - a(n) * (T_now(n-1) - T_now(n))
                end if
        end select

        ! Solve tridiagonal system
        call thomas_solve(a, b, c, d, T_new, n)
    end subroutine solve_implicit

    subroutine thomas_solve(a, b, c, d, x, n)
        integer, intent(in) :: n
        real(8), dimension(n), intent(in) :: a, b, c, d
        real(8), dimension(n), intent(out) :: x
        real(8), dimension(n) :: c_prime, d_prime
        integer :: idx
        real(8) :: m

        c_prime(1) = c(1) / b(1)
        d_prime(1) = d(1) / b(1)

        do idx = 2, n - 1
            m = b(idx) - a(idx) * c_prime(idx-1)
            c_prime(idx) = c(idx) / m
            d_prime(idx) = (d(idx) - a(idx) * d_prime(idx-1)) / m
        end do

        m = b(n) - a(n) * c_prime(n-1)
        d_prime(n) = (d(n) - a(n) * d_prime(n-1)) / m

        x(n) = d_prime(n)
        do idx = n - 1, 1, -1
            x(idx) = d_prime(idx) - c_prime(idx) * x(idx+1)
        end do
    end subroutine thomas_solve

end program transient_fdm_solver