when simulate plasmonic, it's divergent

[ying06]

I want to do some nano photonics plasmonic simulation， the code below is a simple test, it works if I used dielectric in the simuation, but when I use AgUp or Mask material in the simulation, it become divergent， why it happens？ how can I fix it？ The OS is windows 10, version 0.036, with matlab API。
`% Tutorials / radar cross section of a metal sphere

% (C) 2012-2014 Thorsten Liebig thorsten.liebig@uni-due.de

close all
clear
clc

%% setup the simulation
physical_constants;
unit = 1.0e-9; % all length in nm
drawingunit = unit .* 1.0e9;
wlunit = unit;

% size of the simulation box
xmin = -300 .* drawingunit;
xmax = 300 .* drawingunit;
zmin = -600 .* drawingunit;
zmax = 600 .* drawingunit;

Rxmin = xmin - 20.*drawingunit;
Rymin = Rxmin;
Rxmax = xmax + 20.*drawingunit;
Rymax = Rxmax;
Rzmin = zmin - 20.*drawingunit;
Rzmax = zmax + 20.*drawingunit;

mesh_res_xy = 10.0 .* drawingunit;
mesh_res_z = 10.0 .* drawingunit;

%% setup FDTD parameter & excitation function

%FDTD setup
Lam0 = 365 .* wlunit;
Lam_min = 320 .* wlunit;
Lam_max = 410 .* wlunit;
f0 = c0 ./ Lam0; %center frequency
f_min = c0 ./ Lam_max;
f_max = c0 ./ Lam_min;
f_BW = f_max - f_min; %bandwidth
f_start = f_min; % start frequency
f_stop = f_max; % stop frequency

AgUp.eps_R = 1;
AgUp.f0 = 1.85 .* f0; %plasma frequency of the drude material
AgUp.gamma = 0.5 .* f0;
AgUp.relaxTime = 1 ./ AgUp.gamma; %relaxation time (smaller number results in greater losses, set to 0 to disable)

Mask.eps_R = 4.5;
Mask.f0 = 2.0 .* f0; %plasma frequency of the drude material
Mask.gamma = 3.5 .* f0;
Mask.relaxTime = 1 ./ Mask.gamma; %relaxation time (smaller number results in greater losses, set to 0 to disable)

FDTD = InitFDTD('NrTS', 50000, 'EndCriteria', 1e-05);
FDTD = SetGaussExcite(FDTD, f0, f_BW);
BC = [1 1 1 1 1 1] * 3; % set boundary conditions
% BC = [1 1 0 0 2 2];
FDTD = SetBoundaryCond(FDTD, BC);

%% setup CSXCAD geometry & mesh
% max_res = c0 / f_stop / unit / 50; % cell size: lambda/20
CSX = InitCSX();

%create mesh
mesh.x = Rxmin:mesh_res_xy:Rxmax;
mesh.y = Rxmin:mesh_res_xy:Rxmax;
mesh.z = Rzmin:mesh_res_z:Rzmax;

nx = length(mesh.x);
ny = length(mesh.y);
nz = length(mesh.z);

% plane wave excitation
k_dir = [0, 0, 1];

E_dir = [0 1 0];

CSX = AddPlaneWaveExcite(CSX, 'plane_wave', k_dir, E_dir, f0);
% CSX = AddExcitation(CSX, 'plance', 0, [0 1 0]);
start = [xmin xmin zmin];
stop = [xmax xmax zmax];
CSX = AddBox(CSX, 'plane_wave', 0, start, stop);

%% apply AgUp %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
CSX = AddMaterial( CSX, 'RO3010' );
CSX = SetMaterialProperty( CSX, 'RO3010', 'Epsilon', 4.0);
start = [mesh.x(1), mesh.y(1), mesh.z(end)];
stop = [mesh.x(end), mesh.y(end), mesh.z(round(nz/2))];
CSX = AddBox(CSX, 'RO3010', 10, start, stop);

CSX = AddLorentzMaterial(CSX, 'AgUp');
CSX = SetMaterialProperty(CSX, 'AgUp', 'Epsilon', AgUp.eps_R,...
'EpsilonPlasmaFrequency', AgUp.f0, 'EpsilonRelaxTime', AgUp.relaxTime);
% start = [mesh.x(1), mesh.y(1), mesh.z(end-60)];
% stop = [mesh.x(end), mesh.y(end), mesh.z(end-40)];
% CSX = AddBox(CSX, 'AgUp', 10, start, stop);

CSX = AddLorentzMaterial(CSX, 'Mask');
CSX = SetMaterialProperty(CSX, 'Mask', 'Epsilon', Mask.eps_R,...
'EpsilonPlasmaFrequency', Mask.f0, 'EpsilonRelaxTime', Mask.relaxTime);

% start = [mesh.x(round(nx/2) + 7), mesh.y(round(ny/2) + 7), mesh.z(1+60)];
% stop = [mesh.x(round(nx/2) + 17), mesh.y(round(ny/2) + 17), mesh.z(1+80)];
% CSX = AddBox(CSX, 'Mask', 10, start, stop);
%
% start = [mesh.x(1), mesh.y(1), zmin + 120.*drawingunit];
% stop = [mesh.x(end), mesh.y(round(ny/2) - 7), zmin + 160.*drawingunit];
% CSX = AddBox(CSX, 'Mask', 10, start, stop);

%% dump boxes
% CSX = AddDump(CSX, 'Et');
CSX = AddDump(CSX, 'Ef_ra', 'DumpType', 10, 'FileType', 1, 'SubSampling', '2,2,2', 'Frequency', f0);
start = [mesh.x(5), mesh.y(round(ny ./ 2)), mesh.z(10)];
stop = [mesh.x(end - 5), mesh.y(round(ny ./ 2)), mesh.z(end - 10)];
CSX = AddBox(CSX, 'Ef_ra', 0, start, stop);

mesh = AddPML(mesh, 8);

CSX = DefineRectGrid(CSX, unit, mesh);

%% prepare simulation folder
Sim_Path = 'Sim_TFSF';
Sim_CSX = 'Sim_TFSF.xml';

if (exist(Sim_Path, 'dir'))
rmdir(Sim_Path, 's');
end

mkdir(Sim_Path);

rmdir(Sim_Path, 's'); % clear previous directory
[status, message, messageid] = mkdir(Sim_Path); % create empty simulation folder

%% write openEMS compatible xml-file
WriteOpenEMS([Sim_Path '/' Sim_CSX], FDTD, CSX);

%% show the structure
% CSXGeomPlot( [Sim_Path '/' Sim_CSX] );

%% plot the drude type material dependency
% LamList = linspace(Lam_min, Lam_max, 201);
% f = c0 ./ LamList;
% % f = linspace(f_min, f_max, 201);
% w = 2 * pi .* f;
% epsr = AgUp.eps_R .* (1 - (2 * pi * AgUp.f0)^2 ./ (w.^2 - 1j .* w ./ AgUp.relaxTime));
% figure
% plot(LamList, real(epsr), 'Linewidth', 2);
% hold on
% grid on
% plot(LamList, imag(epsr), 'r--', 'Linewidth', 2);
% ylim([-10 AgUp.eps_R])
% % l=legend('\Re \epsilon_r','\Im \epsilon_r','\Re \mue_r','\Im \mue_r');
% l = legend('$\Re{\varepsilon_r}$', '$\Im{\varepsilon_r}$');
% set(l, 'Interpreter', 'latex', 'Fontsize', 12)
%
% LamList = linspace(Lam_min, Lam_max, 201);
% f = c0 ./ LamList;
% % f = linspace(f_min, f_max, 201);
% w = 2 * pi .* f;
% epsr = Mask.eps_R .* (1 - (2 * pi * Mask.f0)^2 ./ (w.^2 - 1j .* w ./ Mask.relaxTime));
% figure
% plot(LamList, real(epsr), 'Linewidth', 2);
% hold on
% grid on
% plot(LamList, imag(epsr), 'r--', 'Linewidth', 2);
% ylim([-10 Mask.eps_R])
% % l=legend('\Re \epsilon_r','\Im \epsilon_r','\Re \mue_r','\Im \mue_r');
% l = legend('$\Re{\varepsilon_r}$', '$\Im{\varepsilon_r}$');
% set(l, 'Interpreter', 'latex', 'Fontsize', 12)
% return
%% run openEMS
RunOpenEMS(Sim_Path, Sim_CSX);

%%
disp('Use Paraview to display the elctric fields dumped by openEMS');

%% plot Field Res
[field, mesh_h5] = ReadHDF5Dump([Sim_Path '/Ef_ra.h5']);

Ey = squeeze(field.FD.values{1}(:, :, :, 2));

X = mesh_h5.lines{1};
Z = mesh_h5.lines{3};

figure
imagesc(X, Z, real(Ey).')
colorbar;
xlabel('x(nm)');
ylabel('z(nm)');
axis xy;
axis image;
MyFormatFinal`

[thliebig]

I think this would be better as a discussion not an issue... because usually thinks like this are just a user error.
Here it looks like you created an overlap between dispersive material and the PML? That does not work. Please check if this is the case.