% FD1D_HEAT_EXPLICIT: Finite difference solution of 1D heat equation.
%
%  Discussion:
%
%    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 explicit forward 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) - 2 U(X,T) + U(X+dx,T) )
%      ------------------  = F(X,T) + k *  ------------------------------------
%               dt                                   dx * dx
%
%    or, assuming we have solved for all values of U at time T, we have
%
%      U(X,T+dt) = U(X,T) + cfl * ( U(X-dx,T) - 2 U(X,T) + U(X+dx,T) ) + dt * F(X,T) 
%
%    Here "cfl" is the Courant-Friedrichs-Loewy coefficient:
%
%      cfl = k * dt / dx / dx
%
%
%  Parameters:
%
%    Input, integer XN, the number of points to use in the spatial dimension.
%
%    Input, integer TN, the number of equally spaced points in the time dimension.
%
%    Output, real H(XN,TN), the computed solution.
%
%    Output, real X(XN), the spatial points.
%
%    Output, real T(TN), the time points.
%
  

function h = boundary_conditions ( xn, x, t, h )
%% BOUNDARY_CONDITIONS evaluates the boundary conditions.
  h(1)  = 90.0;
  h(xn) = 70.0;
end


function h = initial_condition ( xn, x, t )
%% INITIAL_CONDITION evaluates the initial condition.
  h(1:xn) = 50.0;
end


function value = rhs ( x, t )
%% RHS evaluates the right hand side of the differential equation.
  value = x;
  value = 0.0;
end


%  XN is the number of equally spaced nodes to use between 0 and 1.
%
  xn = 21;
%
%  TN is the number of equally spaced time points between 0 and 10.0.
%
  tn = 201;
%
%  Running the code produces an array H of temperatures H(t,x),
%  and vectors x and t.
%


%  Heat coefficient.
%
  k = 0.002;
%
%  X values;
%
  xmin = 0.0;
  xmax = 1.0;
  x = linspace ( xmin, xmax, xn );
  dx = ( xmax - xmin ) / ( xn - 1 );
%
%  T values;
%
  tmin = 0.0;
  tmax = 80.0;
  dt = ( tmax - tmin ) / ( tn - 1 );
  t = linspace ( tmin, tmax, tn );
%
%  Check the CFL condition, have processor 0 print out its value,
%  and quit if it is too large.
%
  cfl = k * dt / dx / dx;

  if ( 0.5 <= cfl )
     error ( 'Error: CFL condition failed. 0.5 <= K * dT / dX / dX' );
  end
%
%  Compute and save initial values.
%
  h = initial_condition ( xn, x, t(1) );
  h = boundary_conditions ( xn, x, t(1), h );

  H = zeros ( tn, xn );
  H(1,1:xn) = h(1:xn);
%
%  Compute the values of H at the next time, based on current data.
%
  L = 1 : xn - 2;
  C = 2 : xn - 1;
  R = 3 : xn;

  for i = 1 : tn - 1

    h(C) = h(C) + cfl * ( h(L) - 2.0 * h(C) + h(R) ) + dt * rhs ( x(C), t(i) );

    h = boundary_conditions ( xn, x, t(i+1), h );

    H(i+1,1:xn) = h(1:xn);

  end


write("position: ",x);

for i = 1 : tn - 1
   p = H(i,1:xn);
   write(i,":  ", p);
end