clear; clc; close all;

%% =========================
% 参数
%% =========================
L = 1.0;
Nx = 256;
dx = L / Nx;

c = 1.0;
CFL = 2.0;   % 👉 故意设大！演示"无条件稳定"
nu = CFL;

dt = CFL * dx / c;

T_end = 0.5;
Nt = round(T_end / dt);

x = linspace(0, L, Nx);

%% =========================
% 初始条件（多频）
%% =========================
u0 = sin(2*pi*5*x) + 0.5*sin(2*pi*20*x) + 0.2*sin(2*pi*60*x);

u = u0;

%% =========================
% 频率
%% =========================
k = (0:Nx-1);
theta = 2*pi*k/Nx;

%% =========================
% 构造BTCS矩阵（周期边界）
%% =========================
e = ones(Nx,1);

A = spdiags([-nu/2*e, e, nu/2*e], [-1 0 1], Nx, Nx);

% 周期边界
A(1,end) = -nu/2;
A(end,1) =  nu/2;

%% =========================
% 图像
%% =========================
figure('Color','w','Position',[100 100 1400 800]);

for n = 1:Nt
    
    %% ===== BTCS推进 =====
    u = A \ u';   % 解线性方程
    u = u';
    
    %% =========================
    % 1️⃣ 物理空间
    %% =========================
    subplot(2,2,1)
    plot(x, u, 'b','LineWidth',1.5);
    title(['Physical Space (BTCS, t = ', num2str(n*dt,'%.2f'), ')']);
    axis([0 1 -2 2]); grid on;
    
    %% =========================
    % 2️⃣ FFT
    %% =========================
    U = fft(u);
    spectrum = abs(U)/Nx;
    
    subplot(2,2,2)
    stem(k(1:Nx/2), spectrum(1:Nx/2),'filled');
    title('Fourier Spectrum');
    xlabel('Wavenumber'); ylabel('Amplitude');
    grid on;
    
    %% =========================
    % 3️⃣ 放大因子 |G|
    %% =========================
    G = 1 ./ (1 + 1i*nu*sin(theta));
    
    subplot(2,2,3)
    plot(theta, abs(G),'r','LineWidth',2); hold on;
    yline(1,'k--','|G|=1');
    hold off;
    
    title(['|G| (BTCS, CFL=',num2str(CFL),')']);
    xlabel('\theta'); ylabel('|G|');
    grid on;
    
    %% =========================
    % 4️⃣ 复平面
    %% =========================
    subplot(2,2,4)
    
    t_circle = linspace(0,2*pi,200);
    plot(cos(t_circle), sin(t_circle),'k--'); hold on;
    
    plot(real(G), imag(G),'r','LineWidth',1.5);
    
    hold off;
    axis equal;
    xlim([-1.5 1.5]); ylim([-1.5 1.5]);
    
    title('Complex Plane (BTCS)');
    xlabel('Re(G)'); ylabel('Im(G)');
    grid on;
    
    drawnow;
end