diff --git a/FDTD/extensions/operator_ext_lorentzmaterial.cpp b/FDTD/extensions/operator_ext_lorentzmaterial.cpp
index 1e6475d..c272204 100644
--- a/FDTD/extensions/operator_ext_lorentzmaterial.cpp
+++ b/FDTD/extensions/operator_ext_lorentzmaterial.cpp
@@ -52,9 +52,10 @@ bool Operator_Ext_LorentzMaterial::BuildExtension()
 	unsigned int numLines[3] = {m_Op->GetNumberOfLines(0),m_Op->GetNumberOfLines(1),m_Op->GetNumberOfLines(2)};
 	CSPropLorentzMaterial* mat = NULL;
 
-	double w_plasma;
+	double w_plasma,t_relax;
 	bool b_pos_on;
 	double L_D[3], C_D[3];
+	double R_D[3], G_D[3];
 	vector<double> v_int[3];
 	vector<double> v_ext[3];
 	vector<double> i_int[3];
@@ -73,6 +74,7 @@ bool Operator_Ext_LorentzMaterial::BuildExtension()
 				for (int n=0; n<3; ++n)
 				{
 					L_D[n]=0;
+					R_D[n]=0;
 					coord[0] = m_Op->GetDiscLine(0,pos[0]);
 					coord[1] = m_Op->GetDiscLine(1,pos[1]);
 					coord[2] = m_Op->GetDiscLine(2,pos[2]);
@@ -82,18 +84,24 @@ bool Operator_Ext_LorentzMaterial::BuildExtension()
 					if ((mat = prop->ToLorentzMaterial()))
 					{
 						w_plasma = mat->GetEpsPlasmaFreqWeighted(n,coord) * 2 * PI;
-						if (w_plasma)
+						if (w_plasma>0)
 						{
 							b_pos_on = true;
 							m_volt_ADE_On = true;
 							L_D[n] = 1/(w_plasma*w_plasma*m_Op->EC_C[n][index]);
 						}
+						t_relax = mat->GetEpsRelaxTimeWeighted(n,coord);
+						if ((t_relax>0) && m_volt_ADE_On)
+						{
+							R_D[n] = L_D[n]/t_relax;
+						}
 					}
 				}
 
 				for (int n=0; n<3; ++n)
 				{
 					C_D[n]=0;
+					G_D[n]=0;
 					coord[0] = m_Op->GetDiscLine(0,pos[0],true);
 					coord[1] = m_Op->GetDiscLine(1,pos[1],true);
 					coord[2] = m_Op->GetDiscLine(2,pos[2],true);
@@ -103,12 +111,17 @@ bool Operator_Ext_LorentzMaterial::BuildExtension()
 					if ((mat = prop->ToLorentzMaterial()))
 					{
 						w_plasma = mat->GetMuePlasmaFreqWeighted(n,coord) * 2 * PI;
-						if (w_plasma)
+						if (w_plasma>0)
 						{
 							b_pos_on = true;
 							m_curr_ADE_On = true;
 							C_D[n] = 1/(w_plasma*w_plasma*m_Op->EC_L[n][index]);
 						}
+						t_relax = mat->GetMueRelaxTimeWeighted(n,coord);
+						if ((t_relax>0) && m_curr_ADE_On)
+						{
+							G_D[n] = C_D[n]/t_relax;
+						}
 					}
 				}
 
@@ -117,16 +130,26 @@ bool Operator_Ext_LorentzMaterial::BuildExtension()
 					for (unsigned int n=0; n<3; ++n)
 					{
 						v_pos[n].push_back(pos[n]);
-						v_int[n].push_back(1);
-						i_int[n].push_back(1);
 						if (L_D[n]>0)
-							v_ext[n].push_back(dT/L_D[n]	*m_Op->GetVI(n,pos[0],pos[1],pos[2]));
+						{
+							v_int[n].push_back((2*L_D[n]-dT*R_D[n])/(2*L_D[n]+dT*R_D[n]));
+							v_ext[n].push_back(dT/(L_D[n]+dT*R_D[n]/2)*m_Op->GetVI(n,pos[0],pos[1],pos[2]));
+						}
 						else
+						{
+							v_int[n].push_back(1);
 							v_ext[n].push_back(0);
+						}
 						if (C_D[n]>0)
-							i_ext[n].push_back(dT/C_D[n]	*m_Op->GetIV(n,pos[0],pos[1],pos[2]));
+						{
+							i_int[n].push_back((2*C_D[n]-dT*G_D[n])/(2*C_D[n]+dT*G_D[n]));
+							i_ext[n].push_back(dT/(C_D[n]+dT*G_D[n]/2)*m_Op->GetIV(n,pos[0],pos[1],pos[2]));
+						}
 						else
+						{
+							i_int[n].push_back(1);
 							i_ext[n].push_back(0);
+						}
 //						cerr << v_int[n].back() << " " << v_ext[n].back() << " " << i_int[n].back() << " " << i_ext[n].back() << endl;
 					}
 				}
diff --git a/matlab/examples/other/Metamaterial_PlaneWave_Drude.m b/matlab/examples/other/Metamaterial_PlaneWave_Drude.m
index 108bf6f..db72b3e 100644
--- a/matlab/examples/other/Metamaterial_PlaneWave_Drude.m
+++ b/matlab/examples/other/Metamaterial_PlaneWave_Drude.m
@@ -1,7 +1,9 @@
 %%%%%%%%%%%%%%%%%%%%%%%
 % example demonstrating double drude meta-material
 %
-% author: Thorsten Liebig @ 2010
+% tested with openEMS v0.0.28
+%
+% author: Thorsten Liebig @ 2010,2012
 %%%%%%%%%%%%%%%%%%%%%%%
 
 close all
@@ -17,27 +19,27 @@ Settings.LogFile = 'openEMS.log';
 pic_size = round([1400 1400/4]); %define the animation picture size
 
 %simulation domain setup (in mm)
-length = 4000;
-width = 100;    
-mesh_res = 5;        % mesh resolution
+length = 500;
+width = 10;
+mesh_res = 0.5;        % mesh resolution
 height = 3*mesh_res; % hight is ony 3 lines with PEC (top/bottom) --> quasi 2D
 
 %FDTD setup
-f0 = 0.5e9;         %center frequency
+f0 = 5e9;         %center frequency
 f_BW = f0/sqrt(2);  %bandwidth
+MTM.eps_R = 1;
+MTM.mue_R = 1;
 MTM.f0 = f0;        %plasma frequency of the drude material
-MTM.length = 1000;  %length of the metamaterial
+MTM.relaxTime = 5e-9; %relaxation time (smaller number results in greater losses, set to 0 to disable)
+MTM.length = 250;  %length of the metamaterial
 N_TS = 5e4;         %number of timesteps
 endCriteria = 1e-5; %stop simulation if signal is at -50dB
 
 %constants
-EPS0 = 8.85418781762e-12;
-MUE0 = 1.256637062e-6;
+physical_constants
 
 %% define openEMS options %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-openEMS_opts = '';
-% openEMS_opts = [openEMS_opts ' --disable-dumps'];
-openEMS_opts = [openEMS_opts ' --engine=fastest'];
+openEMS_opts = '-vvv';
 
 Sim_Path = 'MTM_PW_Drude';
 Sim_CSX = 'MTM_PW_Drude.xml';
@@ -51,8 +53,7 @@ if (postproc_only==0)
 
     %% setup FDTD parameter & excitation function %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     FDTD = InitFDTD(N_TS,endCriteria,'OverSampling',10);
-    FDTD = SetGaussExcite(FDTD,f0,f0/sqrt(2));
-    % FDTD = SetSinusExcite(FDTD,f0*sqrt(2));
+    FDTD = SetGaussExcite(FDTD,0,2*f0);
     BC = [1 1 0 0 2 2];
     FDTD = SetBoundaryCond(FDTD,BC);
 
@@ -70,9 +71,9 @@ if (postproc_only==0)
     CSX = AddBox(CSX,'excite',0 ,start,stop);
 
     %% apply drude material %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-    CSX = AddProperty(CSX,'LorentzMaterial','drude');
-    CSX = SetPropertyArgs(CSX,'LorentzMaterial','drude','PlasmaFrequency','Epsilon',MTM.f0);
-    CSX = SetPropertyArgs(CSX,'LorentzMaterial','drude','PlasmaFrequency','Mue'    ,MTM.f0);
+    CSX = AddLorentzMaterial(CSX,'drude');
+    CSX = SetMaterialProperty(CSX,'drude','Epsilon',MTM.eps_R,'EpsilonPlasmaFrequency',MTM.f0,'EpsilonRelaxTime',MTM.relaxTime);
+    CSX = SetMaterialProperty(CSX,'drude','Mue',MTM.mue_R,'MuePlasmaFrequency',MTM.f0,'MueRelaxTime',MTM.relaxTime);
     start=[mesh.x(1)   mesh.y(1)   -MTM.length/2];
     stop =[mesh.x(end) mesh.y(end)  MTM.length/2];
     CSX = AddBox(CSX,'drude', 10 ,start,stop);
@@ -91,6 +92,22 @@ if (postproc_only==0)
 
 end
 
+%% plot the drude type material dependency
+f = linspace(0.1*f0,2*f0,501);
+w = 2*pi*f;
+epsr = MTM.eps_R * (1 - (2*pi*MTM.f0)^2./( w.^2 - 1j*w./MTM.relaxTime ));
+muer = MTM.mue_R * (1 - (2*pi*MTM.f0)^2./( w.^2 - 1j*w./MTM.relaxTime ));
+plot(f,real(epsr),'Linewidth',2);
+hold on
+grid on
+plot(f,imag(epsr),'r--','Linewidth',2);
+plot(f,real(muer),'c-.','Linewidth',2);
+plot(f,imag(muer),'m-.','Linewidth',2);
+ylim([-10 MTM.eps_R])
+% l=legend('\Re \epsilon_r','\Im \epsilon_r','\Re \mue_r','\Im \mue_r');
+l=legend('$\Re\{\varepsilon_r\}$','$\Im\{\varepsilon_r\}$','$\Re\{\mu_r\}$','$\Im\{\mu_r\}$');
+set(l,'Interpreter','latex','Fontsize',12)
+
 %% plot E-fields
 freq = [f0/sqrt(2) f0 f0*sqrt(2)];
 field = ReadHDF5FieldData([Sim_Path '/Et.h5']);