まずは基底クラスですよね！有界領域で空間1Dなので, こんな感じかな. ほとんどpSolver1DBase.hpp と同じだが

// bSolver1DBase.hpp    Solve some equation using some method. Bounded region
#pragma once
#include <iostream>
#include <filesystem>
#include <cmath>
#include <vector>
#include <fstream>
#include "../include/pvector.hpp"
template <typename T>
class bSolver1DBase
{
protected:
   double C = 1.0, K=1.0; // Cは速度、Kは拡散係数
   double CFL, DN; // CFL条件と拡散数
   double dt, dx, time, xmax = 1.0;
   int nt = 0; // 時間ステップのカウンタ
   double LBC_value = 0.0; // 下端の境界条件の値
   double RBC_value = 0.0; // 上端の境界条件の値
   class Point
   {
   public:
       Point *pre = nullptr; // ポインタを追加して前のポイントへのリンクを作成
       Point *next = nullptr;
       double X;
       pvector<T> U;     // U_n
       Point() : U(3) {} // デフォルトコンストラクタでサイズ3のpvectorを初期化
       auto &operator[](int index) { return U[index]; }
       const auto &operator[](int index) const { return U[index]; }
   };
   std::vector<Point> Data;
   std::ofstream ofs;
public:
   bSolver1DBase(const size_t &nx, const double &dt, const std::function<double(double)> &U0) : dt(dt)
   {
       time = 0.0;         // 初期時間
       dx = xmax / nx;     // 空間刻み幅の計算
       CFL = C * dt / dx;  // CFL条件の計算
       DN = C * dt /dx/dx; // 拡散数の計算
       Data.resize(nx+1);  // 空間分割数に基づいて初期化
       for (size_t i = 0; auto &P : Data)
       {
           P.pre = &Data[i - 1];
           P.next = &Data[i + 1];
           P[nt] = P[nt - 1] = P[nt - 2] = U0(P.X = i * dx); // 初期条件の適用
           i++;
       }
       std::cout << "Booting bSolver1DBase with nx=" << Data.size() << ", dt=" << dt << std::endl;
   }
   std::string what() const
   {
       return "nSolver1DBase: nx=" + std::to_string(Data.size()) + ", dt=" + std::to_string(dt) +
              ", CFL=" + std::to_string(CFL) + " , DN=" + std::to_string(DN)  + ", time=" + std::to_string(time);
   }
   double& LBC() { return LBC_value; } // 下端の境界条件
   double& RBC() { return RBC_value; } // 上端の境界条件
   bool Write(const std::string &filename)
   {
       if (ofs) ofs.close(); // 既存のファイルストリームを閉じる
       std::filesystem::path datafile = getenv("HOME");
       datafile /= filename;
       std::filesystem::create_directories(datafile.parent_path());
       ofs.open(datafile);
       if (!ofs) throw std::ios_base::failure("Failed to open file: " + datafile.string());
       return Write();
   }
   bool Write()
   {
       try
       {
           ofs << std::endl
               << "# nt= " << nt << " time= " << time << " CFL= " << CFL << " NX= " << Data.size() << std::endl;
           for (auto &v : Data)
               ofs << v.X << " " << v[nt] << std::endl;
       }
       catch (const std::exception &e)
       {
           std::cerr << "Error writing to file: " << e.what() << std::endl;
           return false;
       }
       std::cout << "time[" << time << "]written to file successfully." << std::endl;
       return true;
   }
   virtual void Step(double U_W = 0.0) = 0;
   virtual void Initialize(void *parm = nullptr) {};
   double get_time() const { return time; }
};

まあ名前が違うところ以外は, ここいらで拡散係数などを定義．ここはフツーのstd::vectorで, でも格子点が0,1,...,NX+1 であると, 区間幅は0,1,...,NX になるから,　ここの数が1個違うだけ． 人によっては，pSolver1DBase と bSolver1DBase の共通のご先祖様  Solver1DBase から派生した方がいいと思うかもだな．

あ，ここは！先ほど別のプログラムで作成した pvector.hpp ですね．こういう，よく使う汎用プログラムは，共通のフォルダーを作って，そこに入れておくのが良いです．

include/ というフォルダーが多いかな？ lib/ もありだが. 

で, この辺は，有界領域の問題なので, 上下で境界条件を設定したくなるかもなあ．というものです．このあたりは, まあ今は使いませんが，のちのち, コンストラクタとStep()の間に作業が入るようになるかもしれないので，念の為に作ってあります．

基底クラスがよくできてれば，派生クラスは簡単です

// heat1D_explicit.hpp
#pragma once
#include "bSolver1DBase.hpp"
class Heat1DExplicit : public bSolver1DBase<double>
{
public:
   Heat1DExplicit(const size_t &nx, const double &dt, const std::function<double(double)> &U0 = [](double x)
                  { return 0.; }) : bSolver1DBase<double>(nx, dt, U0)
   {
       std::cout << "Starting Heat1DExplicit solver with nx=" << Data.size() << ", dt=" << dt << std::endl;
   }
   void Step() override
   {
       Data[0][nt + 1] = LBC_value; // 下端の境界条件
       Data[Data.size() - 1][nt + 1] = RBC_value; // 上端の境界条件
       for (size_t i = 1; i < Data.size() - 1; ++i)
           Data[i][nt + 1] = (1.-2.*DN)* Data[i][nt] + DN * (Data[i - 1][nt] + Data[i + 1][nt]);
       time += dt;
       nt++;
   }
};
// heat1D_explicit.hpp

(2)式を書いただけ，ですね．これでmain()は

// main.cpp
#include "heat1D_explicit.hpp"
int main()
{
   size_t nx = 20;          // 空間分割数 100(text), 40, 20, 8
   double dt = 0.1 / nx /nx;    // 時間刻み幅
   size_t ntmax = 5. / dt;  // 最大時間ステップ数
   size_t ntsave = .005 / dt; // データ保存間隔
   std::string filename = "Data/test04/CDE_20.dat";
   Heat1DExplicit mySolver(nx, dt);
   mySolver.LBC() = 1.0; // 下端の境界条件
   mySolver.Initialize();
   mySolver.Write(filename);
   for (size_t nt = 0; nt < ntmax; nt++)
   {
       mySolver.Step(); // 時間ステップを実行
       if ((nt + 1) % ntsave == 0)
           mySolver.Write();
   }
   return 0;
}

ここは(3)にならないように書きました！境界条件を左が1に設定！細かめにデータ保存したら, 左みたいになっちまったorz (一応収束していることは，わかる．でも，こんな図を課長に見せたら「最後どうなってんのか見えねえぞ，帰れ」となってしまう)ので，ここいらへんを改造して右図みたいにしましたとさ. 

改造はこんな感じにしてみた:

#include <set>
   ...
   std::set<int> tset;
   tset.insert(round(0.005/dt));tset.insert(round(0.010/dt));tset.insert(round(0.020/dt));
   tset.insert(round(0.050/dt));tset.insert(round(0.10/dt));tset.insert(round(0.20/dt));
   tset.insert(round(0.50/dt));tset.insert(round(1/dt));tset.insert(round(2/dt));
   tset.insert(round(5/dt));
   ...
   for (size_t nt = 0; nt < ntmax; nt++)
       ...
       if (tset.contains(nt+1))
           mySolver.Write();

なんで整数にしてるんだ? という点ですが, 最初double でやったら，2とかが1.999999999999999 になって, 「2.00000000000000001じゃないから保存しないもん」になったからです．このように, double 数値を A=B とか A>B とかで比較すると，たまにハマる．というのは鉄則なので, 覚えておくと良いです．

陽解法は1Dだと解けますね！でもnt=20000までかかっています.

Eigen/SPARSEの使い方をimplicitスキームの構築で学びましょう. 　　(注意: EigenのVersion 3.4.0を利用してくださいね. EigenのVersionの確認方法は, 猿に聞け)

// heat1D_implicit.hpp
#pragma once
#include "bSolver1DBase.hpp"
#include <Eigen/Sparse>
class Heat1DImplicit : public bSolver1DBase<double>
{
   Eigen::SparseLU<Eigen::SparseMatrix<double>> solver;
public:
   Heat1DImplicit(const size_t &nx, const double &dt, const std::function<double(double)> &U0 = [](double x)
                                                      { return 0.; }) : bSolver1DBase<double>(nx, dt, U0)
   {
       std::cout << "Starting Heat1DImplicit solver with nx=" << Data.size() << ", dt=" << dt << std::endl;
   }
   void Initialize(void *parm = nullptr) override
   {
       Eigen::SparseMatrix<double> coef; // スパース行列を使用して係数行列を表現
       coef.resize(Data.size(), Data.size());
       for (int i = 0; i < Data.size(); i++)
       {
           if ((i==0)||(i==Data.size()-1)) // 境界条件の処理
               coef.insert(i,i) = 1;
           else
           {
               coef.insert(i, i - 1) = -(DN + 0.5 * CFL); // 左隣の成分
               coef.insert(i, i) = 1. + 2. * DN;          // 対角成分だにゃ
               coef.insert(i, i + 1) = -(DN - 0.5 * CFL); // 右隣の成分
           }
       }
       coef.makeCompressed(); // スパース行列を圧縮       
       solver.compute(coef);
       if (solver.info()!=Eigen::Success)
              throw std::runtime_error("Error::Heat1DImplicit::Initialize::Can not solve fuck.");
       std::cout << "Heat1DImplicit solver initialized with nx=" << Data.size() << ", dt=" << dt << std::endl;
   }
   // 時間ステップを実行するメソッド
   void Step() override
   {
       Eigen::VectorXd b(Data.size());
       for (int i = 0; i < Data.size(); i++)
           if (i == 0)                     b[i] = LBC_value; // 下端の境界条件
           else if (i == Data.size() - 1)  b[i] = RBC_value; // 上端の境界条件
           else                            b[i] = Data[i][nt]; // 現在の値を設定
       auto result=solver.solve(b); // 連立方程式を解く
       time = (++nt) * dt;
       for (int i = 0; i < Data.size(); i++) Data[i][nt] = result[i];
   }
};
// heat1D_implicit.hpp

これでEigen/Sparseを読み込みます. まずはLU分解を入れておくところを用意. で, Initialize()で，(4)の行列を作成します．まずはこうやってSparse行列と, そのサイズを設定します.

ここでは1行目とN行目の対角成分を1に設定. この見慣れない｜｜は, 二つの条件の論理OR（どちらかがtrueならtrue）を意味しています. なお, AND(どっちもtrueならtrue)は && です．

なんで | と & ではないのか？試しに書いてみると動くんだが？という向きのために説明しておくと, | は bit-OR  つまり二進数で 0101 と 1001 なら　1101 になる，という計算になってしまいます．&も bit-ANDで, 0101 & 1001 = 0001 になります．で，たまに論理OR, 論理ANDと結果が一致していると，動いている気がする，というわけです．この説明がわからん，という人は，黙って ||   と && と書いてくださいね！

この辺では, i行目の成分として, i-1列目, i列目, i+1列目を設定します．Eigen/Sparseは, こうやって「値を設定」したあと, このmakeCompressedで計算機メモリーを節約します．さて，設定した行列をLU分解しておけば計算が速くなります．ここで必ずやっておきましょう！で, LU分解はいつでも可能ではなく，失敗する場合がありますので，こうやって失敗を検出しておきましょう．

さて, 時間ステップ Step() では, 計算を行います．
では連立方程式\(\mathsf{A}\boldsymbol{x}=\boldsymbol{b}\) \(\mathbf{b}\)ベクトルを設定します.で，これで解けるので, 回収しておきましょう．

ここらへんで解いている最中には，計算機の使用メモリーが上昇します．どれくらい上昇するかは行列に依存しますが，まあ数倍くらいで済むことが多い感じです．

　陰解法なので, 時間ステップをいくらでも大きくとれる. この例ではnt=20でt=5まで到達．

こんな感じ:

#!/bin/bash
#   fig.sh: 偏微分方程式II: 拡散方程式
. draw.sh
. standards.sh
COLOR=1
X_AXIS='\fIX\sb\fR1\eb'
X_RANGE='0 1 0.25'
Y_AXIS='\fIu'
Y_RANGE='-0.1 1.4 0.5'
OPT="--pen-colors=1=DarkRed:2=MediumSeaGreen:3=Blue:4=Magenta:5=Khaki --symbol-line-width 0.0020"
#(c,κ)の解析解のデータが，ぽっちい〜
echo -e "1 1\n1 1e-1\n1 1e-2\n"|gawk '{
NX=100; dx=1/NX; omega=$1/$2;
for(i=0;i<NX+1;i++)
{
   X=i*dx;
   Y=(exp(omega*(X-1))-1)/(exp(-omega)-1);# exp(omega)なんて書いたらオーバーフローするき！
   printf("%f %f\n",X,Y);
}
printf("\n");
}'  >   $TMP/steady
#凡例
cat > $TMP/legend.dat <<@
#m=Solid:Black,S=Omit
0.05 1.3    l:implicit finite difference
0.05 1.15   l:\fIc\fR\r8=\r81, \fI\*D\sbt\eb\fR\r8=\r80.25

0.25    0.825       lb:\*k\r8=\r81

0.75    0.95        rt:\*k\r8=\r80.1

0.98    0.75
0.8     0.7     r:\*k\r8=\r810\sp-2\ep

#m=Dot:Black,S=Omit
0.01    0.5
0.10    0.5         l:\fIN\sbx\eb\fR\r8=\r820, CFL\r8=\r85

#m=Solid:Green,S=Omit
0.01    0.375
0.10    0.375       l:\fIN\sbx\eb\fR\r8=\r840, CFL\r8=\r810

#m=Solid:Magenta,S=Omit
0.01    0.25
0.10    0.25        l:\fIN\sbx\eb\fR\r8=\r880, CFL\r8=\r820
@
cat>>$TMP/steady<<@
0.01    0.125
0.10    0.125       l:theory
@
LEGEND=$TMP/legend.dat
LEGEND_OPT="--line-width 0.003"
cat>$TMP/nls<<@
$TMP/steady             Solid:Cyan:0.007        Omit
ADi_20_K1.dat.last      Dot:Black:0.003     Omit
ADi_20_K1e-1.dat.last   Dot:Black:0.003     Omit
ADi_20_K1e-2.dat.last   Dot:Black:0.003     Omit
ADi_40_K1e-2.dat.last   Solid:Green:0.006   Omit
ADi_80_K1e-2.dat.last   Solid:Magenta:0.003 Omit
@
NAMAE=(     `cat $TMP/nls|grep -v ^#|gawk '{printf("%s ",$1);}'` )
LINE=(      `cat $TMP/nls|grep -v ^#|gawk '{printf("%s ",$2);}'` )
SYMBOL=(    `cat $TMP/nls|grep -v ^#|gawk '{printf("%s ",$3);}'` )
draw ADi.pdf ${NAMAE[@]}

計算プログラムよりも，こっちがむしろ大変だなwww