To accurately evaluate a central difference derivative (which is $O(\Delta x^2)$) exactly at a boundary node ($i=0$), we mathematically imagine a fictitious "ghost cell" ($i=-1$) located outside the physical domain. For a standard Neumann condition like $\frac{\partial T}{\partial x} = 0$, we set the central difference: $\frac{T_1 - T_{-1}}{2\Delta x} = 0 \implies T_{-1} = T_1$ This algebraic ghost value $T_{-1}$ is then substituted directly into the main governing PDE evaluated at the boundary node $i=0$. This crucial technique allows the highly accurate central difference scheme to be maintained throughout the entire domain without dropping to a less accurate, asymmetric forward/backward difference ($O(\Delta x)$) at the edges.

For a square computational grid ($\Delta x = \Delta y$), the continuous PDE algebraically simplifies to a basic averaging formula: $T_{i,j} = \frac{T_{i+1,j} + T_{i-1,j} + T_{i,j+1} + T_{i,j-1}}{4}$ Applying this formula to every unknown internal node yields a large, sparse system of simultaneous linear algebraic equations ($[A]\{T\} = \{b\}$). This system can be solved directly (e.g., Gaussian elimination) or iteratively (e.g., via the Gauss-Seidel method, historically termed Liebmann's method for this specific application).

MOL discretizes only the spatial variables (using standard finite differences or finite elements) while leaving time $t$ continuous. This powerful technique algebraically reduces the single continuous PDE to a massive, coupled system of Ordinary Differential Equations (ODEs) exclusively with respect to time. These resulting ODEs can then be passed to highly optimized, off-the-shelf adaptive ODE solvers (like Runge-Kutta-Fehlberg or Gear's methods for stiff systems) which automatically handle time-stepping stability. It provides an elegant, rigorous separation of spatial discretization errors from time integration errors.

The fundamental mathematical stability of any linear finite difference scheme is rigorously evaluated using Von Neumann stability analysis. By substituting a discrete Fourier series component $T_j^n = \xi^n e^{ik(j\Delta x)}$ into the difference equation, we determine the amplification factor $\xi$. For the numerical solution to remain bounded (stable) and not violently oscillate to infinity, the strict condition $|\xi| \le 1$ must be satisfied for all possible frequencies $k$. This analysis mathematically proves why explicit schemes have severe time-step restrictions, while implicit schemes are often unconditionally stable regardless of $\Delta t$.

For the 1D heat equation, the explicit finite difference approximation yields an algebraic stepping formula based on the dimensionless parameter $r = k\frac{\Delta t}{(\Delta x)^2}$. The Von Neumann mathematical stability criterion requires that the coefficient of $T_i^n$ remains non-negative, rigorously implying: $r \le \frac{1}{2} \implies \Delta t \le \frac{(\Delta x)^2}{2k}$ If this strict constraint is violated by even a fraction, the numerical solution will immediately become violently unstable, producing infinite non-physical oscillations.

The Crank-Nicolson method is a widely used implicit numerical scheme that averages the spatial difference approximations equally at the current known time step and the unknown future time step: $\frac{T_i^{n+1} - T_i^n}{\Delta t} = k \left[ \frac{1}{2} \left( \frac{T_{i+1}^{n} - 2T_i^{n} + T_{i-1}^{n}}{(\Delta x)^2} \right) + \frac{1}{2} \left( \frac{T_{i+1}^{n+1} - 2T_i^{n+1} + T_{i-1}^{n+1}}{(\Delta x)^2} \right) \right]$ This balanced averaging makes the method unconditionally stable (via Von Neumann analysis) and highly second-order accurate in both space and time ($O(\Delta x^2, \Delta t^2)$). This mathematical property allows for much larger time steps than explicit methods without any risk of numerical instability or explosion.

For 2D transient parabolic equations, simultaneously solving a fully implicit large 2D sparse matrix system at every single time step is computationally crushing. The ADI method elegantly splits the physical time step $\Delta t$ in half. In the first mathematical half-step, it solves implicitly in the x-direction and explicitly in the y-direction. In the second half-step, it precisely reverses the directions. This brilliant formulation yields simple 1D tridiagonal matrices for every row/column that can be solved extraordinarily rapidly using the Thomas algorithm, while maintaining unconditional stability.

While discrete numerical methods are common for complex domains, the fundamental 1D wave equation has a famous exact analytical solution formulated by d'Alembert. It mathematically proves that any arbitrary solution can be expressed precisely as the superposition of two rigid waves traveling in opposite directions at constant speed $c$ without changing shape: $u(x,t) = F(x - ct) + G(x + ct)$ Where $F$ represents the arbitrary wave profile traveling strictly to the right, and $G$ represents the wave traveling strictly to the left. The specific geometric forms of $F$ and $G$ are determined entirely by the given initial position and velocity conditions at $t=0$.

The CFL condition dictates the absolute maximum allowable numerical time step $\Delta t$ relative to the spatial grid size $\Delta x$ for explicit hyperbolic solvers. It is universally expressed in terms of the dimensionless Courant number $C$: $C = \frac{c\Delta t}{\Delta x}$ Where $c$ is the physical wave speed. For numerical stability, it mathematically requires $C \leq 1$. Physically, this means the numerical domain of dependence (the information traveling across the discrete grid at speed $\Delta x / \Delta t$) must encompass the true physical domain of dependence (the physical wave traveling at speed $c$). If the physical wave outruns the numerical grid calculation, the simulation catastrophically fails.