写一段FDTD的matlab代码

以下是一个简单的二维FDTD Matlab代码，用于计算一个矩形波导中的电磁场分布：

% FDTD simulation of a rectangular waveguide

% Define simulation parameters
dx = 0.01;  % spatial step size
dt = 0.001; % time step size
T = 1;      % total simulation time
Lx = 1;     % waveguide length
Ly = 0.5;   % waveguide width

% Define material properties
mu0 = 4*pi*10^-7;      % permeability of free space
epsilon0 = 8.854e-12;  % permittivity of free space
epsilonr = 2.5;        % relative permittivity of waveguide
mur = 1;               % relative permeability of waveguide

% Calculate material constants
c = 1/sqrt(mu0*epsilon0);    % speed of light in free space
epsilon = epsilonr*epsilon0; % permittivity of waveguide
mu = mur*mu0;                % permeability of waveguide

% Calculate grid dimensions
Nx = round(Lx/dx);
Ny = round(Ly/dx);

% Initialize electric and magnetic fields
Ex = zeros(Nx,Ny);
Ey = zeros(Nx,Ny);
Hz = zeros(Nx,Ny);

% Define source function
f = @(t) sin(2*pi*10*t);

% Main simulation loop
for n = 1:round(T/dt)
    
    % Update electric fields
    for i = 1:Nx
        for j = 1:Ny-1
            Ex(i,j) = Ex(i,j) + (dt/(dx*epsilon))*...
                (Hz(i,j) - Hz(i,j+1));
        end
    end
    for i = 1:Nx-1
        for j = 1:Ny
            Ey(i,j) = Ey(i,j) - (dt/(dx*epsilon))*...
                (Hz(i,j) - Hz(i+1,j));
        end
    end
    
    % Apply source function
    Ey(Nx/2,Ny/2) = Ey(Nx/2,Ny/2) + f(n*dt);
    
    % Update magnetic field
    for i = 1:Nx-1
        for j = 1:Ny-1
            Hz(i,j) = Hz(i,j) + (dt/(dx*mu))*...
                (Ex(i,j) - Ex(i,j+1) + Ey(i+1,j) - Ey(i,j));
        end
    end
    
    % Plot electric field
    imagesc(Ex');
    colormap(jet);
    colorbar;
    drawnow;
end

此代码使用了基本的FDTD算法，通过在每个网格点处计算电场和磁场分量，并使用Maxwell方程组进行递推计算，得到电磁场的时变分布。在这个例子中，我们在波导中心的一个网格点处施加了一个10 Hz的正弦波形源，然后计算了电场的分布情况。