请写出三维有限时域差分的matlab代码

以下是一个简单的三维有限时域差分法的 MATLAB 代码示例：

% 定义计算参数
dx = 0.01; % 空间分辨率
dt = dx / c0; % 时间分辨率，c0为光速常数
N = 100; % 空间步长
T = 1000; % 时间步长

% 定义空间和时间网格
Ex = zeros(N, N, N); % 电场x分量
Ey = zeros(N, N, N); % 电场y分量
Ez = zeros(N, N, N); % 电场z分量
Hx = zeros(N, N, N); % 磁场x分量
Hy = zeros(N, N, N); % 磁场y分量
Hz = zeros(N, N, N); % 磁场z分量

% 定义电磁波源函数
source = zeros(T, 1);
source(1) = 1;

% 迭代计算电磁场变化
for t = 1:T
    % 更新电场
    Ex(2:N-1, 2:N-1, 2:N-1) = Ex(2:N-1, 2:N-1, 2:N-1) + dt / eps0 / dx * (Hz(2:N-1, 2:N-1, 2:N-1) - Hz(2:N-1, 1:N-2, 2:N-1) - Hy(2:N-1, 2:N-1, 2:N-1) + Hy(2:N-1, 2:N-1, 1:N-2));
    Ey(2:N-1, 2:N-1, 2:N-1) = Ey(2:N-1, 2:N-1, 2:N-1) + dt / eps0 / dx * (Hx(2:N-1, 2:N-1, 2:N-1) - Hx(2:N-1, 2:N-1, 1:N-2) - Hz(2:N-1, 2:N-1, 2:N-1) + Hz(1:N-2, 2:N-1, 2:N-1));
    Ez(2:N-1, 2:N-1, 2:N-1) = Ez(2:N-1, 2:N-1, 2:N-1) + dt / eps0 / dx * (Hy(2:N-1, 2:N-1, 2:N-1) - Hy(1:N-2, 2:N-1, 2:N-1) - Hx(2:N-1, 2:N-1, 2:N-1) + Hx(2:N-1, 1:N-2, 2:N-1));
    % 更新磁场
    Hx(2:N-1, 2:N-1, 2:N-1) = Hx(2:N-1, 2:N-1, 2:N-1) - dt / mu0 / dx * (Ez(2:N-1, 2:N-1, 2:N-1) - Ez(2:N-1, 1:N-2, 2:N-1) - Ey(2:N-1, 2:N-1, 2:N-1) + Ey(2:N-1, 2:N-1, 1:N-2));
    Hy(2:N-1, 2:N-1, 2:N-1) = Hy(2:N-1, 2:N-1, 2:N-1) - dt / mu0 / dx * (Ex(2:N-1, 2:N-1, 2:N-1) - Ex(2:N-1, 2:N-1, 1:N-2) - Ez(2:N-1, 2:N-1, 2:N-1) + Ez(1:N-2, 2:N-1, 2:N-1));
    Hz(2:N-1, 2:N-1, 2:N-1) = Hz(2:N-1, 2:N-1, 2:N-1) - dt / mu0 / dx * (Ey(2:N-1, 2:N-1, 2:N-1) - Ey(1:N-2, 2:N-1, 2:N-1) - Ex(2:N-1, 2:N-1, 2:N-1) + Ex(2:N-1, 1:N-2, 2:N-1));
    % 添加源函数
    Ex(N/2, N/2, N/2) = Ex(N/2, N/2, N/2) + source(t);
    % 绘制电磁波传播动画
    slice(Ex, [], [], N/2);
    shading interp;
    axis equal;
    axis tight;
    drawnow;
end

该代码实现了一个简单的三维电磁波传播动画。在每个时间步长中，通过更新电磁场和磁场的分量，可以计算出电磁场在空间和时间上的变化情况。同时，添加源函数可以模拟电磁波的发射和接收过程。