#! /usr/bin/env python
#
def fd1d_heat_implicit ( a, x_num, x, t, dt, cfl, rhs_fun, bc_fun, u ):

#*****************************************************************************80
#
## FD1D_HEAT_IMPLICIT: Finite difference solution of 1D heat equation.
#
#  Discussion:
#
#    FD1D_HEAT_IMPLICIT solves the 1D heat equation with an implicit method.
#
#    This program solves
#
#      dUdT - k * d2UdX2 = F(X,T)
#
#    over the interval [A,B] with boundary conditions
#
#      U(A,T) = UA(T),
#      U(B,T) = UB(T),
#
#    over the time interval [T0,T1] with initial conditions
#
#      U(X,T0) = U0(X)
#
#    The code uses the finite difference method to approximate the
#    second derivative in space, and an implicit backward Euler approximation
#    to the first derivative in time.
#
#    The finite difference form can be written as
#
#      U(X,T+dt) - U(X,T)                  ( U(X-dx,T+dt) - 2 U(X,T+dt) + U(X+dx,T+dt) )
#      ------------------ = F(X,T+dt) + k *  --------------------------------------
#               dt                                   dx * dx
#
#    so that we have the following linear system for the values of U at time T+dt:
#
#            -     k * dt / dx / dx   * U(X-dt,T+dt)
#      + ( 1 + 2 * k * dt / dx / dx ) * U(X,   T+dt)
#            -     k * dt / dx / dx   * U(X+dt,T+dt)
#      =               dt             * F(X,   T+dt)
#      +                                U(X,   T)
#
#  Licensing:
#
#    This code is distributed under the GNU LGPL license.
#
#  Modified:
#
#    14 April 2017
#
#  Author:
#
#    John Burkardt
#
#  Parameters:
#
#    Input, real A(X_NUM,X_NUM), the system matrix.
#
#    Input, integer X_NUM, the number of nodes.
#
#    Input, real X(X_NUM), the node coordinates.
#
#    Input, real T, the current time.
#
#    Input, real DT, the size of the time step.
#
#    Input, real CFL, the Courant-Friedrichs-Loewy coefficient.
#
#    Input, f = RHS_FUN ( x_num, x, t ), returns in F the right hand side
#    forcing function at every non-boundary node.
#
#    Input, hbc = BC_FUN ( x_num, x, t, h ), returns in HBC a copy of the
#    input solution H, after imposing Dirichlet boundary conditions.
#
#    Input, real U(X_NUM), the solution values at the old time.
#
#    Output, real U(X_NUM), the solution values at the new time.
#
  import numpy as np
#
#  Compute b, the right hand side of the system.
#
  fvec = rhs_fun ( x_num, x, t )

  b = u.copy ( )
  for i in range ( 1, x_num - 1 ):
    b[i] = b[i] + dt * fvec[i]
#
#  Solve A*u=b.
#
  u = np.linalg.solve ( a, b )
#
#  Impose boundary conditions on U.
#
  u = bc_fun ( x_num, x, t, u )

  return u