clear; clc; close all;

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

alpha = 0.01;
dt = 0.4 * dx^2 / alpha;

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

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

%% =========================================================
% 初始条件选择
%% =========================================================
IC_type =  4;

switch IC_type
    
    case 1
        u = exp(-100*(x-0.5).^2);
        title_str = '高斯（解析解可用）';
        
    case 2
        u = double(x>0.3 & x<0.7);
        title_str = '方波（无解析解）';
        
    case 3
        u = sin(pi*x);
        title_str = '正弦（解析解可用）';
        
    case 4
        rng(1);
        u = 2*rand(size(x)) - 0.5;
        title_str = '随机噪声（无解析解）';
        
end

u0 = u;

% 👉 存储 x-t 数据
U = zeros(Nt+1, Nx);
U(1,:) = u;

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

% ===== 左图：1D =====
subplot(1,2,1)
h_num = plot(x,u,'b-','LineWidth',2); hold on;
h_exact = plot(x,u,'r--','LineWidth',1.5);

legend('数值解','解析解');
xlabel('x'); ylabel('u');
grid on;
axis([0 L -1 1.5]);

htitle = title(title_str);

% ===== 右图：x-t =====
subplot(1,2,2)
h_img = imagesc(x, t_all, U);
set(gca,'YDir','normal');
xlabel('x'); ylabel('t');
title('x-t 演化图（扩散过程）');
colorbar;

%% ---------------- 时间推进 ----------------
for n = 1:Nt
    
    % 数值解
    u(2:end-1) = u(2:end-1) + alpha*dt/dx^2 * ...
        (u(3:end) - 2*u(2:end-1) + u(1:end-2));
    
    % 边界
    u(1) = 0;
    u(end) = 0;
    
    t = n*dt;
    
    % 👉 存储
    U(n+1,:) = u;
    
    %% ===== 精确解 =====
    switch IC_type
        
        case 1
            u_exact = 1/sqrt(1+400*alpha*t) * ...
                exp(-100*(x-0.5).^2/(1+400*alpha*t));
            
        case 3
            u_exact = exp(-alpha*pi^2*t)*sin(pi*x);
            
        otherwise
            u_exact = nan(size(x));
    end
    
    %% ===== 刷新 =====
    if mod(n,2)==0
        
        % 左图
        subplot(1,2,1)
        set(h_num,'YData',u);
        set(h_exact,'YData',u_exact);
        htitle.String = [title_str,'， t = ', num2str(t,'%.3f')];
        
        % 右图
        subplot(1,2,2)
        set(h_img,'CData',U);
        
        drawnow limitrate
        pause(0.001)
    end
end