Eddy current field solving in PDE tool

Hi!

A band of long ideal conductors placed in a magnetic field which is changing sinusoidal in time. In the conductors evolve of line dipole character eddy currents. It can be described with a complex permeability. The task is to model the system and describe how the complex permeability depends on frequency.

I have to change frequency between 10kHz and 10MHz. I can't make it in Matlab, because matlab handles the coeff's of the PDE parameters as strings, so it seems thjat I dont use the freq variable.

My code:

model:

% This script is written and read by pdetool and should NOT be edited.
% There are two recommended alternatives:
% 1) Export the required variables from pdetool and create a MATLAB script
%    to perform operations on these.
% 2) Define the problem completely using a MATLAB script. See
%    https://www.mathworks.com/help/pde/examples.html for examples
%    of this approach.
function pdemodel(freq,B_0)
[pde_fig,ax]=pdeinit;
pdetool('appl_cb',1);
pdetool('snapon','on');
set(ax,'DataAspectRatio',[1 1 1]);
set(ax,'PlotBoxAspectRatio',[1.6434567901234571 1 1666.6666666666665]);
set(ax,'XLimMode','auto');
set(ax,'YLim',[-0.00040000000000000002 0.001]);
set(ax,'XTickMode','auto');
set(ax,'YTickMode','auto');
pdetool('gridon','on');


% Geometry description:
pdeellip(0.000275000,0.0004763,0.00025000,0.00025000,...
0,'E2');
pdeellip(0,0,0.00025000,0.00025000,...
0,'E1');
pderect([0 0.000275000 0 0.0004763],'R1');
set(findobj(get(pde_fig,'Children'),'Tag','PDEEval'),'String','R1')

% Boundary conditions:
pdetool('changemode',0)
pdesetbd(8,...
'dir',...
1,...
'1',...
'0')
pdesetbd(7,...
'dir',...
1,...
'1',...
'0')
pdesetbd(6,...
'dir',...
1,...
'1',...
'0')
pdesetbd(5,...
'dir',...
1,...
'1',...
'0')
pdesetbd(4,...
'dir',...
1,...
'1',...
'0')
pdesetbd(3,...
'dir',...
1,...
'1',...
'0')
pdesetbd(2,...
'dir',...
1,...
'1',...
'0')
pdesetbd(1,...
'dir',...
1,...
'1',...
'0')

% Mesh generation:
setappdata(pde_fig,'Hgrad',1.3);
setappdata(pde_fig,'refinemethod','regular');
setappdata(pde_fig,'jiggle',char('on','mean',''));
setappdata(pde_fig,'MesherVersion','preR2013a');
pdetool('initmesh')
pdetool('refine')
pdetool('refine')

% PDE coefficients:
pdeseteq(1,...
'1./(x*pi*4e-7)!1./(x*pi*4e-7)!1./(x*pi*4e-7)',...
'0!-2i*58.82*pi*freq!-2i*58.82*pi*freq',...
'0!-B_0*x*58.82*2*pi*freq!-B_0*x*58.82*2*pi*freq',...
'1.0!1.0!1.0',...
'0:10',...
'0.0',...
'0.0',...
'[0 100]')
setappdata(pde_fig,'currparam',...
['1./(x*pi*4e-7)!1./(x*pi*4e-7)!1./(x*pi*4e-7)   ';...
'0!-2i*58.82*pi*freq!-2i*58.82*pi*freq          ';...
'0!-B_0*x*58.82*2*pi*freq!-B_0*x*58.82*2*pi*freq';...
'1.0!1.0!1.0                                    '])

% Solve parameters:
setappdata(pde_fig,'solveparam',...
char('0','41256','10','pdeadworst',...
'0.5','longest','0','1E-4','','fixed','Inf'))

% Plotflags and user data strings:
setappdata(pde_fig,'plotflags',[1 1 1 1 1 1 1 1 0 0 0 1 1 0 0 0 0 1]);
setappdata(pde_fig,'colstring','');
setappdata(pde_fig,'arrowstring','');
setappdata(pde_fig,'deformstring','');
setappdata(pde_fig,'heightstring','');

Script which varies the frequency

%Logarithmic scale of frequencies
B_0=100;
i=3;
while(i<=6)
i=i+0.25;
freq=10^(1+i);
pdemodel(freq,B_0); 
end
    

pde_fig=findobj(allchild(0),'flat','Tag','PDETool');
u = get(findobj(pde_fig,'Tag','PDEPlotMenu'),'UserData');
h=findobj(get(pde_fig,'Children'),'flat','Tag','PDEMeshMenu');
hp=findobj(get(h,'Children'),'flat','Tag','PDEInitMesh');
he=findobj(get(h,'Children'),'flat','Tag','PDERefine');
ht=findobj(get(h,'Children'),'flat','Tag','PDEMeshParam');
p=get(hp,'UserData');
e=get(he,'UserData');
t=get(ht,'UserData');
hbound = findobj(get(pde_fig,'Children'),'flat','Tag','PDEBoundMenu');
g = get(hbound, 'userdata');

Can sy please help me?

Thanks

Vitya