clear; clc; close all;

%% ---------------- 参数 ----------------
L = 1.0;
Nx = 200;
dx = L / Nx;

c = 1.0;
CFL = 0.1;                
dt = CFL * dx / c;

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

x = linspace(0, L, Nx);
t = (0:Nt)*dt;

%% ---------------- 初始条件 ----------------
u0 = exp(-100 * (x - 0.3).^2);

%% ---------------- 精确解 ----------------
u_exact_fun = @(x,t) exp(-100 * (mod(x - c*t, L) - 0.3).^2);

%% ---------------- 初始化 ----------------
u_ftcs = u0;
u_up   = u0;
u_lax  = u0;
u_btcs = u0;

%% ---------------- BTCS矩阵 ----------------
r = c * dt / dx;
e = ones(Nx,1);

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

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

%% ---------------- 存储 ----------------
U_ftcs = zeros(Nt+1, Nx);
U_up   = zeros(Nt+1, Nx);
U_lax  = zeros(Nt+1, Nx);
U_btcs = zeros(Nt+1, Nx);

U_ftcs(1,:) = u_ftcs;
U_up(1,:)   = u_up;
U_lax(1,:)  = u_lax;
U_btcs(1,:) = u_btcs;

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

% ===== 左图：解对比 =====
subplot(2,3,[1 4])
h1 = plot(x, u_ftcs, 'r-', 'LineWidth',1.5); hold on;
h2 = plot(x, u_up,   'b-', 'LineWidth',1.5);
h3 = plot(x, u_lax,  'g-', 'LineWidth',1.5);
h4 = plot(x, u_btcs, 'k-', 'LineWidth',1.5);

% 精确解
u_exact = u_exact_fun(x,0);
h5 = plot(x, u_exact, 'm--', 'LineWidth',2);

legend('FTCS（不稳定）','Upwind','Lax','BTCS','Exact');

xlabel('x'); ylabel('u');
grid on;
axis([0 L -0.5 1.5]);

htitle = title('初始化');

% ===== 右侧四个 x-t 图 =====
subplot(2,3,2)
contourf(x, t, U_ftcs, 20, 'LineColor','none');
title('FTCS');
xlabel('x'); ylabel('t'); colorbar;

subplot(2,3,3)
contourf(x, t, U_up, 20, 'LineColor','none');
title('Upwind');
xlabel('x'); ylabel('t'); colorbar;

subplot(2,3,5)
contourf(x, t, U_lax, 20, 'LineColor','none');
title('Lax');
xlabel('x'); ylabel('t'); colorbar;

subplot(2,3,6)
contourf(x, t, U_btcs, 20, 'LineColor','none');
title('BTCS');
xlabel('x'); ylabel('t'); colorbar;

%% ---------------- 时间推进 ----------------
for n = 1:Nt
    
    %% ===== FTCS（不稳定）=====
    u_new = u_ftcs;
    u_new(2:Nx-1) = u_ftcs(2:Nx-1) ...
        - CFL/2 * (u_ftcs(3:Nx) - u_ftcs(1:Nx-2));
    
    u_new(1)  = u_new(Nx-1);
    u_new(Nx) = u_new(2);
    u_ftcs = u_new;
    
    %% ===== Upwind =====
    u_new = u_up;
    u_new(2:Nx) = u_up(2:Nx) ...
        - CFL * (u_up(2:Nx) - u_up(1:Nx-1));
    u_new(1) = u_new(Nx);
    u_up = u_new;
    
    %% ===== Lax =====
    u_new = u_lax;
    u_new(2:Nx-1) = 0.5*(u_lax(3:Nx) + u_lax(1:Nx-2)) ...
        - CFL/2*(u_lax(3:Nx) - u_lax(1:Nx-2));
    
    u_new(1)  = u_new(Nx-1);
    u_new(Nx) = u_new(2);
    u_lax = u_new;
    
    %% ===== BTCS（隐式）=====
    u_btcs = A \ u_btcs';
    u_btcs = u_btcs';
    
    %% ===== 存储 =====
    U_ftcs(n+1,:) = u_ftcs;
    U_up(n+1,:)   = u_up;
    U_lax(n+1,:)  = u_lax;
    U_btcs(n+1,:) = u_btcs;
    
    %% ===== 可视化 =====
    if mod(n,5)==0
        
        % 左图更新
        subplot(2,3,[1 4])
        set(h1,'YData',u_ftcs);
        set(h2,'YData',u_up);
        set(h3,'YData',u_lax);
        set(h4,'YData',u_btcs);
        
        % 精确解更新
        u_exact = u_exact_fun(x, n*dt);
        set(h5,'YData', u_exact);
        
        htitle.String = ['t = ', num2str(n*dt,'%.2f'), ...
                         ' | CFL = ', num2str(CFL)];
        
        % 右图更新
        subplot(2,3,2)
        contourf(x, t, U_ftcs, 20, 'LineColor','none');
        
        subplot(2,3,3)
        contourf(x, t, U_up, 20, 'LineColor','none');
        
        subplot(2,3,5)
        contourf(x, t, U_lax, 20, 'LineColor','none');
        
        subplot(2,3,6)
        contourf(x, t, U_btcs, 20, 'LineColor','none');
        
        drawnow limitrate
    end
end