% COMSOL Multiphysics Model M-file
% Generated by COMSOL 3.4 (COMSOL 3.4.0.248, $Date: 2007/10/10 16:07:51 $)

flclear fem

% COMSOL version
clear vrsn
vrsn.name = 'COMSOL 3.4';
vrsn.ext = '';
vrsn.major = 0;
vrsn.build = 248;
vrsn.rcs = '$Name:  $';
vrsn.date = '$Date: 2007/10/10 16:07:51 $';
fem.version = vrsn;

% Geometry
g1=rect2(2,1,'base','corner','pos',[-1,0]);

% Analyzed geometry
clear s
s.objs={g1};
s.name={'R1'};
s.tags={'g1'};

fem.draw=struct('s',s);
fem.geom=geomcsg(fem);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = 'q0';
bnd.ind = [1,1,1,1];
appl.bnd = bnd;
clear equ
equ.C = '1.61*10000/(20*10^(-12))';
equ.init = 300;
equ.k = '386*10^10';
equ.Q = '4.8*10^16*(T2-T)+2.58*10^26*exp(-1*(z*10^4/15.9+1108*t^2))';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = 'q0';
bnd.ind = [1,1,1,1];
appl.bnd = bnd;
clear equ
equ.C = '3.4*10^6/(20*10^(-12))';
equ.init = 300;
equ.k = 0;
equ.Q = '4.8*10^16*(T-T2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.1:1], ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=1    Surface: Temperature [K]', ...
         'axis',[-1.8253012048192772,1.394578313253012,-0.42168674698795205,1.578313253012048]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=1    Surface: Temperature [K]', ...
         'axis',[-1.8253012048192772,1.394578313253012,-0.48356190720328135,1.6401884132273774]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=1    Surface: Temperature [K]', ...
         'axis',[-1.8253012048192772,1.394578313253012,-0.48356190720328135,1.6401884132273774]);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5);

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5);

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5);

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = 'q0';
bnd.ind = [1,1,1,1];
appl.bnd = bnd;
clear equ
equ.C = '1.61*10000/(20*10^(-12))';
equ.init = 300;
equ.k = '386*10^10';
equ.Q = '4.8*10^16*(T2-T)+2.58*10^26*exp(-1*(z*10^4/15.9+1108*t^2))';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = 'q0';
bnd.ind = [1,1,1,1];
appl.bnd = bnd;
clear equ
equ.C = '3.4*10^6/(20*10^(-12))';
equ.init = 300;
equ.k = 0;
equ.Q = '4.8*10^16*(T-T2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.1:1], ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.696751    Surface: Temperature [K]', ...
         'axis',[-1.8253012048192772,1.394578313253012,-0.42168674698795205,1.578313253012048]);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = 'q0';
bnd.ind = [1,1,1,1];
appl.bnd = bnd;
clear equ
equ.C = '1.61*10000/(20*10^(-12))';
equ.init = 300;
equ.k = '386*10^10';
equ.Q = '4.8*10^16*(T2-T)+2.58*10^26*exp(-1*(z*10^4/15.9+1108*t^2))';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = 'q0';
bnd.ind = [1,1,1,1];
appl.bnd = bnd;
clear equ
equ.C = '3.4*10^6/(20*10^(-12))';
equ.init = 300;
equ.k = 0;
equ.Q = '4.8*10^16*(T-T2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.1:1], ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=1    Surface: Temperature [K]', ...
         'axis',[-1.8253012048192772,1.394578313253012,-0.42168674698795205,1.578313253012048]);

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = 'q0';
bnd.ind = [1,1,1,1];
appl.bnd = bnd;
clear equ
equ.C = '1.61*10000/(20*10^(-12))';
equ.init = 300;
equ.k = '386*10^10';
equ.Q = '4.8*10^16*(T2-T)+2.58*10^26*exp(-1*(z*10^4/15.9+1108*t^2))';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = 'q0';
bnd.ind = [1,1,1,1];
appl.bnd = bnd;
clear equ
equ.C = '3.4*10^6/(20*10^(-12))';
equ.init = 300;
equ.k = 0;
equ.Q = '4.8*10^16*(T-T2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.1:1], ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.696722    Surface: Temperature [K]', ...
         'axis',[-1.8253012048192772,1.394578313253012,-0.42168674698795205,1.578313253012048]);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = 'q0';
bnd.ind = [1,1,1,1];
appl.bnd = bnd;
clear equ
equ.C = '1.61*10000/(20*10^(-12))';
equ.init = 300;
equ.k = '386*10^10';
equ.Q = '4.8*10^16*(T2-T)+2.58*10^26*exp(-1*(z*10^4/15.9+1108*t^2))';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = 'q0';
bnd.ind = [1,1,1,1];
appl.bnd = bnd;
clear equ
equ.C = '3.4*10^6/(20*10^(-12))';
equ.init = 300;
equ.k = 0;
equ.Q = '4.8*10^16*(T-T2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:0.11], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.11    Surface: Temperature [K]', ...
         'axis',[-1.8253012048192772,1.394578313253012,-0.42168674698795205,1.578313253012048]);

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = 'q0';
bnd.ind = [1,1,1,1];
appl.bnd = bnd;
clear equ
equ.C = '1.61*10000/(20*10^(-12))';
equ.init = 300;
equ.k = '386*10^10';
equ.Q = '4.8*10^16*(T2-T)+2.58*10^26*exp(-1*(z*10^4/15.9+1108*t^2))';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = 'q0';
bnd.ind = [1,1,1,1];
appl.bnd = bnd;
clear equ
equ.C = '3.4*10^6/(20*10^(-12))';
equ.init = 300;
equ.k = 0;
equ.Q = '4.8*10^16*(T-T2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:0.11], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.11    Surface: Temperature [K]', ...
         'axis',[-1.8253012048192772,1.394578313253012,-0.42168674698795205,1.578313253012048]);

% Animate solution
postmovie(fem, ...
          'tridata',{'T2','cont','internal','unit','K'}, ...
          'trimap','jet(1024)', ...
          'axis',[-1.8253012048192772,1.394578313253012,-0.49005655650123453,1.6466830625253306], ...
          'fps',10);

% Animate solution
postmovie(fem, ...
          'tridata',{'T2','cont','internal','unit','K'}, ...
          'trimap','jet(1024)', ...
          'axis',[-1.8253012048192772,1.394578313253012,-0.49005655650123453,1.6466830625253306], ...
          'fps',10);

% Geometry
g2=rect2(2.8,1,'base','corner','pos',[-1.8,0]);
g3=rect2(3.4,1,'base','corner','pos',[-1.8,0]);
g4=rect2(3.4,0.6,'base','corner','pos',[-1.8,0]);
g4=move(g4,[0.10000000000000009,0.09999999999999998]);
g4=move(g4,[-0.2,0]);
g5=rect2(3.8,0.6,'base','corner','pos',[-1.9,0.1]);
g6=rect2(4,0.6,'base','corner','pos',[-1.9,0.1]);
g6=move(g6,[-0.10000000000000009,0]);

% Analyzed geometry
clear s
s.objs={g6};
s.name={'R1'};
s.tags={'g6'};

fem.draw=struct('s',s);
fem.geom=geomcsg(fem);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5);

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = 'q0';
bnd.ind = [1,1,1,1];
appl.bnd = bnd;
clear equ
equ.C = '1.61*10000/(20*10^(-12))';
equ.init = 300;
equ.k = '386*10^10';
equ.Q = '4.8*10^16*(T2-T)+2.58*10^26*exp(-1*(z*10^4/15.9+1108*t^2))';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = 'q0';
bnd.ind = [1,1,1,1];
appl.bnd = bnd;
clear equ
equ.C = '3.4*10^6/(20*10^(-12))';
equ.init = 300;
equ.k = 0;
equ.Q = '4.8*10^16*(T-T2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:0.11], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.11    Surface: Temperature [K]', ...
         'axis',[-1.180722891566267,2.0391566265060237,-0.4216867469879525,1.5783132530120485]);

% Geometry
g6=move(g6,[0,-0.2]);
g6=move(g6,[0,-0.4]);
g6=move(g6,[0,-0.2]);
g7=rect2('4','0.6','base','corner','pos',{'-2','0'},'rot','0');
g8=rect2('4','0.6','base','corner','pos',{'-2','-0.6'},'rot','0');

% Analyzed geometry
clear s
s.objs={g8};
s.name={'R1'};
s.tags={'g8'};

fem.draw=struct('s',s);
fem.geom=geomcsg(fem);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5);

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5);

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = 'q0';
bnd.ind = [1,1,1,1];
appl.bnd = bnd;
clear equ
equ.C = '1.61*10000/(20*10^(-12))';
equ.init = 300;
equ.k = '386*10^10';
equ.Q = '4.8*10^16*(T2-T)+2.58*10^26*exp(-1*(z*10^4/15.9+1108*t^2))';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = 'q0';
bnd.ind = [1,1,1,1];
appl.bnd = bnd;
clear equ
equ.C = '3.4*10^6/(20*10^(-12))';
equ.init = 300;
equ.k = 0;
equ.Q = '4.8*10^16*(T-T2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:0.11], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.11    Surface: Temperature [K]', ...
         'axis',[-1.1807228915662678,2.0391566265060246,-0.6626506024096395,1.337349397590362]);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'q0','T'};
bnd.T0 = {0,300};
bnd.ind = [2,2,1,2];
appl.bnd = bnd;
clear equ
equ.C = '1.61*10000/(20*10^(-12))';
equ.init = 300;
equ.k = '386*10^10';
equ.Q = '4.8*10^16*(T2-T)+2.58*10^26*exp(-1*(z*10^4/15.9+1108*t^2))';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'q0','T'};
bnd.T0 = {0,300};
bnd.ind = [2,2,1,2];
appl.bnd = bnd;
clear equ
equ.C = '3.4*10^6/(20*10^(-12))';
equ.init = 300;
equ.k = 0;
equ.Q = '4.8*10^16*(T-T2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:0.11], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.11    Surface: Temperature [K]', ...
         'axis',[-1.1807228915662678,2.0391566265060246,-0.6626506024096395,1.337349397590362]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'triz',{'T','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.11    Surface: Temperature [K]   Height: Temperature [K]', ...
         'grid','on');

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'triz',{'T','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.11    Surface: Temperature [K]   Height: Temperature [K]', ...
         'grid','on');

% Plot solution
postplot(fem, ...
         'contdata',{'T','cont','internal','unit','K'}, ...
         'contlevels',20, ...
         'contlabel','off', ...
         'contmap','cool(1024)', ...
         'solnum','end', ...
         'title','Time=0.11    Contour: Temperature [K]', ...
         'axis',[-2,2,-0.6000000238418579,0]);

% Plot solution
postplot(fem, ...
         'arrowdata',{'fluxr_ht','fluxz_ht'}, ...
         'arrowxspacing',15, ...
         'arrowyspacing',15, ...
         'arrowtype','arrow', ...
         'arrowstyle','proportional', ...
         'arrowcolor',[1.0,0.0,0.0], ...
         'solnum','end', ...
         'title','Time=0.11    Arrow: Heat flux', ...
         'axis',[-1.180722891566267,2.0391566265060237,-0.7375949953660804,1.412293790546803]);

% Plot solution
postplot(fem, ...
         'flowdata',{'fluxr_ht','fluxz_ht'}, ...
         'flowcolor',[1.0,0.0,0.0], ...
         'flowlines',20, ...
         'solnum','end', ...
         'title','Time=0.11    Streamline: Heat flux', ...
         'axis',[-1.21287222623212,2.0713059611718765,-0.6626506024096395,1.337349397590362]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.11    Surface: Temperature [K]', ...
         'axis',[-1.21287222623212,2.0713059611718765,-0.6626506024096395,1.337349397590362]);

% Animate solution
postmovie(fem, ...
          'tridata',{'T2','cont','internal','unit','K'}, ...
          'trimap','jet(1024)', ...
          'axis',[-1.180722891566267,2.0391566265060237,-0.7375949953660804,1.412293790546803], ...
          'fps',10);
% COMSOL Multiphysics Model M-file
% Generated by COMSOL 3.4 (COMSOL 3.4.0.248, $Date: 2007/10/10 16:07:51 $)

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'q0','T'};
bnd.T0 = {0,300};
bnd.ind = [2,2,1,2];
appl.bnd = bnd;
clear equ
equ.C = '1.61*10000/(20*10^(-12))';
equ.init = 300;
equ.k = '386*10^10';
equ.Q = '4.8*10^16*(T2-T)+2.58*10^26*exp(-1*(-1*z*10^4/15.9+1108*t^2))';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'q0','T'};
bnd.T0 = {0,300};
bnd.ind = [2,2,1,2];
appl.bnd = bnd;
clear equ
equ.C = '3.4*10^6/(20*10^(-12))';
equ.init = 300;
equ.k = 0;
equ.Q = '4.8*10^16*(T-T2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:0.11], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.11    Surface: Temperature [K]', ...
         'axis',[-1.325335291301092,2.1507446209272865,-0.7045705900525185,1.412293790546803]);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'q0','T'};
bnd.T0 = {0,300};
bnd.ind = [2,2,1,2];
appl.bnd = bnd;
clear equ
equ.C = '1.61*10000/(20*10^(-12))';
equ.init = 300;
equ.k = '386*10^10';
equ.Q = '4.8*10^16*(T2-T)+2.58*10^26*exp(-1*(-1*z*10^4/15.9+1108*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'q0','T'};
bnd.T0 = {0,300};
bnd.ind = [2,2,1,2];
appl.bnd = bnd;
clear equ
equ.C = '3.4*10^6/(20*10^(-12))';
equ.init = 300;
equ.k = 0;
equ.Q = '4.8*10^16*(T-T2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:0.11], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.11    Surface: Temperature [K]', ...
         'axis',[-1.325335291301092,2.1507446209272865,-0.7045705900525185,1.412293790546803]);

% Geometry
g1=rect2('1','0.2','base','corner','pos',{'-2','-0.6'},'rot','0');
g2=rect2('1','0.2','base','corner','pos',{'-0.5','-0.2'},'rot','0');

% Analyzed geometry
clear s
s.objs={g2};
s.name={'R1'};
s.tags={'g2'};

fem.draw=struct('s',s);
fem.geom=geomcsg(fem);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '1.61*10000/(20*10^(-12))';
equ.init = 300;
equ.k = '386*10^10';
equ.Q = '4.8*10^16*(T2-T)+2.58*10^26*exp(-1*(-1*z*10^4/15.9+1108*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '3.4*10^6/(20*10^(-12))';
equ.init = 300;
equ.k = 0;
equ.Q = '4.8*10^16*(T-T2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:0.11], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.11    Surface: Temperature [K]', ...
         'axis',[-1.325335291301092,2.1507446209272865,-0.7045705900525185,1.412293790546803]);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5);

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '1.61*10000/(20*10^(-12))';
equ.init = 300;
equ.k = '386*10^10';
equ.Q = '4.8*10^16*(T2-T)+2.58*10^26*exp(-1*(-1*z*10^4/15.9+1108*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '3.4*10^6/(20*10^(-12))';
equ.init = 300;
equ.k = 0;
equ.Q = '4.8*10^16*(T-T2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:0.11], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.11    Surface: Temperature [K]', ...
         'axis',[-1.2913074005940373,2.1167167302202317,-0.7045705900525185,1.412293790546803]);
% COMSOL Multiphysics Model M-file
% Generated by COMSOL 3.4 (COMSOL 3.4.0.248, $Date: 2007/10/10 16:07:51 $)

% Geometry
g1=rect2('1','0.3','base','corner','pos',{'-0.5','-0.3'},'rot','0');

% Analyzed geometry
clear s
s.objs={g1};
s.name={'R1'};
s.tags={'g1'};

fem.draw=struct('s',s);
fem.geom=geomcsg(fem);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5);

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5);

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '1.61*10000/(20*10^(-12))';
equ.init = 300;
equ.k = '386*10^12';
equ.Q = '4.8*10^16*(T2-T)+2.58*10^26*exp(-1*(-1*z*10^3/15.9+1108*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '3.4*10^6/(20*10^(-12))';
equ.init = 300;
equ.k = 0;
equ.Q = '4.8*10^16*(T-T2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:0.11], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.11    Surface: Temperature [K]', ...
         'axis',[-0.35015774549905915,0.5005670751252549,-0.2607464473277731,0.26846964782205757]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'triz',{'T','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.11    Surface: Temperature [K]   Height: Temperature [K]', ...
         'grid','on');

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.11    Surface: Temperature [K]', ...
         'axis',[-0.5,0.5,-0.30000001192092896,0]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'triz',{'T','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.11    Surface: Temperature [K]   Height: Temperature [K]', ...
         'grid','on');

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.11    Surface: Temperature [K]', ...
         'axis',[-0.5,0.5,-0.30000001192092896,0]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.11    Surface: Temperature [K]', ...
         'axis',[-0.3259752782843548,0.4763846079105506,-0.2611529113025428,0.2688761117968273]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.11    Surface: Temperature [K]', ...
         'axis',[-0.3259752782843548,0.4763846079105506,-0.2611529113025428,0.2688761117968273]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.11    Surface: Temperature [K]', ...
         'axis',[-0.3259752782843548,0.4763846079105506,-0.2611529113025428,0.2688761117968273]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'triz',{'T','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.11    Surface: Temperature [K]   Height: Temperature [K]', ...
         'grid','on');

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.11    Surface: Temperature [K]', ...
         'axis',[-0.5,0.5,-0.30000001192092896,0]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'tridlim',[300.0 3000], ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.11    Surface: Temperature [K]', ...
         'axis',[-0.3259752782843548,0.4763846079105506,-0.2611529113025428,0.2688761117968273]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'tridlim',[300.0 3000], ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.11    Surface: Temperature [K]', ...
         'axis',[-0.3318534206434053,0.4822627502696011,-0.2607464473277731,0.26846964782205757]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'tridlim',[300.0 3000], ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.11    Surface: Temperature [K]', ...
         'axis',[-0.7389115060999085,0.8893208357261043,-0.5253544949026884,0.533077695396973]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'tridlim',[300.0 3000], ...
         'triz',{'T','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.11    Surface: Temperature [K]   Height: Temperature [K]', ...
         'grid','on');

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'triz',{'T','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.11    Surface: Temperature [K]   Height: Temperature [K]', ...
         'grid','on');

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'triz',{'T','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.11    Surface: Temperature [K]   Height: Temperature [K]', ...
         'grid','on');

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.11    Surface: Temperature [K]', ...
         'axis',[-0.5,0.5,-0.30000001192092896,0]);

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.11    Surface: Temperature [K]', ...
         'axis',[-0.7389115060999085,0.8893208357261043,-0.5339334792212267,0.5416566797155113]);

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.11    Surface: Temperature [K]', ...
         'axis',[-0.7389115060999085,0.8893208357261043,-0.5306839621247411,0.5384071626190258]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.11    Surface: Temperature [K]', ...
         'axis',[-0.7389115060999085,0.8893208357261043,-0.5306839621247411,0.5384071626190258]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.11    Surface: Temperature [K]', ...
         'axis',[-0.7389115060999085,0.8893208357261043,-0.5339334792212267,0.5416566797155113]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum',2, ...
         'title','Time=0.01    Surface: Temperature [K]', ...
         'axis',[-0.7389115060999085,0.8893208357261043,-0.5339334792212267,0.5416566797155113]);

% Animate solution
postmovie(fem, ...
          'tridata',{'T2','cont','internal','unit','K'}, ...
          'trimap','jet(1024)', ...
          'axis',[-0.7389115060999085,0.8893208357261044,-0.5339334792212267,0.5416566797155113], ...
          'fps',10);
% COMSOL Multiphysics Model M-file
% Generated by COMSOL 3.4 (COMSOL 3.4.0.248, $Date: 2007/10/10 16:07:51 $)

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5);

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5);

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '1.61*10000/(20*10^(-12))';
equ.init = 300;
equ.k = '386*10^12';
equ.Q = '4.8*10^16*(T2-T)+2.58*10^26*exp(-1*(-1*z*10^3/15.9+1108*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '3.4*10^6/(20*10^(-12))';
equ.init = 300;
equ.k = 0;
equ.Q = '4.8*10^16*(T-T2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.1:2], ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.5840488759893655,0.5840488759893655,-0.5133232075815385,0.2133231956606095]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum',3, ...
         'title','Time=0.2    Surface: Temperature [K]', ...
         'axis',[-0.5500000000000002,0.5500000000000002,-0.5133232075815385,0.2133231956606095]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum',3, ...
         'title','Time=0.2    Surface: Temperature [K]', ...
         'axis',[-0.5500000000000002,0.5500000000000002,-0.5133232075815385,0.2133231956606095]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum',5, ...
         'title','Time=0.4    Surface: Temperature [K]', ...
         'axis',[-0.5500000000000002,0.5500000000000002,-0.5133232075815385,0.2133231956606095]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum',5, ...
         'title','Time=0.4    Surface: Temperature [K]', ...
         'axis',[-0.5500000000000002,0.5500000000000002,-0.5133232075815385,0.2133231956606095]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum',7, ...
         'title','Time=0.6    Surface: Temperature [K]', ...
         'axis',[-0.5500000000000002,0.5500000000000002,-0.5133232075815385,0.2133231956606095]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum',7, ...
         'title','Time=0.6    Surface: Temperature [K]', ...
         'axis',[-0.5500000000000002,0.5500000000000002,-0.5133232075815385,0.2133231956606095]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum',11, ...
         'title','Time=1    Surface: Temperature [K]', ...
         'axis',[-0.5500000000000002,0.5500000000000002,-0.5133232075815385,0.2133231956606095]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum',11, ...
         'title','Time=1    Surface: Temperature [K]', ...
         'axis',[-0.5500000000000002,0.5500000000000002,-0.5133232075815385,0.2133231956606095]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum',16, ...
         'title','Time=1.5    Surface: Temperature [K]', ...
         'axis',[-0.5500000000000002,0.5500000000000002,-0.5133232075815385,0.2133231956606095]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum',16, ...
         'title','Time=1.5    Surface: Temperature [K]', ...
         'axis',[-0.5500000000000002,0.5500000000000002,-0.5133232075815385,0.2133231956606095]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.5500000000000002,0.5500000000000002,-0.5133232075815385,0.2133231956606095]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.5500000000000002,0.5500000000000002,-0.5133232075815385,0.2133231956606095]);

% Animate solution
postmovie(fem, ...
          'tridata',{'T2','cont','internal','unit','K'}, ...
          'trimap','jet(1024)', ...
          'axis',[-0.5500000000000002,0.5499999999999999,-0.5133232075815385,0.2133231956606095], ...
          'fps',10);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '1.61*10000/(20*10^(-12))';
equ.init = 300;
equ.k = '386*10^12';
equ.Q = '4.8*10^16*(T2-T)+2.58*10^26*exp(-1*(-1*z*10^3/15.9+1108*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '3.4*10^6/(20*10^(-12))';
equ.init = 300;
equ.k = 0;
equ.Q = '4.8*10^16*(T-T2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.001:0.2], ...
                'atol',{'0.000010'}, ...
                'rtol',0.0001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.2    Surface: Temperature [K]', ...
         'axis',[-0.5500000000000002,0.5499999999999999,-0.5133232075815385,0.2133231956606095]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum',2, ...
         'title','Time=0.001    Surface: Temperature [K]', ...
         'axis',[-0.5533434650455928,0.5533434650455928,-0.5133232075815385,0.2133231956606095]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum',2, ...
         'title','Time=0.001    Surface: Temperature [K]', ...
         'axis',[-0.5533434650455928,0.5533434650455928,-0.5155318592631255,0.21553184734219655]);
% COMSOL Multiphysics Model M-file
% Generated by COMSOL 3.4 (COMSOL 3.4.0.248, $Date: 2007/10/10 16:07:51 $)

% Animate solution
postmovie(fem, ...
          'tridata',{'T2','cont','internal','unit','K'}, ...
          'trimap','jet(1024)', ...
          'axis',[-0.5533434650455928,0.5533434650455928,-0.5155318592631256,0.21553184734219663], ...
          'solnum',[1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129,130,131,132,133,134,135,136,137,138,139,140,141,142,143,144,145,146,147,148,149,150,151,152,153,154,155,156,157,158,159,160,161,162,163,164,165,166,167,168,169,170,171,172,173,174,175,176,177,178,179,180,181,182,183,184,185,186,187,188,189,190,191,192,193,194,195,196,197,198,199,200,201], ...
          'fps',10);

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum',2, ...
         'title','Time=0.001    Surface: Temperature [K]', ...
         'axis',[-0.5533434650455928,0.5533434650455928,-0.5155318592631256,0.21553184734219663]);

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum',2, ...
         'title','Time=0.001    Surface: Temperature [K]', ...
         'axis',[-0.5533434650455928,0.5533434650455928,-0.5155318592631256,0.21553184734219663]);

% Animate solution
postmovie(fem, ...
          'tridata',{'T','cont','internal','unit','K'}, ...
          'trimap','jet(1024)', ...
          'axis',[-0.5533434650455928,0.5533434650455928,-0.5155318592631256,0.21553184734219663], ...
          'fps',10);
% COMSOL Multiphysics Model M-file
% Generated by COMSOL 3.4 (COMSOL 3.4.0.248, $Date: 2007/10/10 16:07:51 $)

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '1.61*10000/(20*10^(-12))';
equ.init = 300;
equ.k = '386*10^12';
equ.Q = '4.8*10^16*(T2-T)+2.58*10^26*exp(-1*(-1*z*10^3/15.9+1108*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '3.4*10^6/(20*10^(-12))';
equ.init = 300;
equ.k = 0;
equ.Q = '4.8*10^16*(T-T2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.001:0.2], ...
                'atol',{'0.000010'}, ...
                'rtol',0.0001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0 (2)    Surface: Temperature [K]', ...
         'axis',[-0.5875993250835259,0.5875993250835259,-0.5155318592631256,0.21553184734219663]);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '1.61*10000/(20*10^(-12))';
equ.init = 300;
equ.k = '386*10^12';
equ.Q = '4.8*10^16*(T2-T)+2.58*10^26*exp(-1*(-1*z*10^3/15.9+1108*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '3.4*10^6/(20*10^(-12))';
equ.init = 300;
equ.k = 0;
equ.Q = '4.8*10^16*(T-T2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.0001:0.2], ...
                'atol',{'0.0000010'}, ...
                'rtol',0.00001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0 (2)    Surface: Temperature [K]', ...
         'axis',[-0.5656773619361735,0.5656773619361735,-0.5155318592631256,0.21553184734219663]);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5);

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '1.61*10000/(20*10^(-12))';
equ.init = 300;
equ.k = '386*10^12';
equ.Q = '4.8*10^16*(T2-T)+2.58*10^26*exp(-1*(-1*z*10^3/15.9+1108*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '3.4*10^6/(20*10^(-12))';
equ.init = 300;
equ.k = 0;
equ.Q = '4.8*10^16*(T-T2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.001:0.2], ...
                'atol',{'0.000010'}, ...
                'rtol',0.0001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0 (2)    Surface: Temperature [K]', ...
         'axis',[-0.5875993250835259,0.5875993250835259,-0.5155318592631256,0.21553184734219663]);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5);

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '1.61*10000/(20*10^(-12))';
equ.init = 300;
equ.k = '386*10^12';
equ.Q = '4.8*10^16*(T2-T)+2.58*10^26*exp(-1*(-1*z*10^3/15.9+1108*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '3.4*10^6/(20*10^(-12))';
equ.init = 300;
equ.k = 0;
equ.Q = '4.8*10^16*(T-T2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.001:0.2], ...
                'atol',{'0.000010'}, ...
                'rtol',0.0001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.2    Surface: Temperature [K]', ...
         'axis',[-0.5875993250835259,0.5875993250835259,-0.5155318592631256,0.21553184734219663]);

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '1.61*10000/(20*10^(-12))';
equ.init = 300;
equ.k = '386*10^12';
equ.Q = '4.8*10^16*(T2-T)+2.58*10^26*exp(-1*(-1*z*10^3/15.9+1108*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '3.4*10^6/(20*10^(-12))';
equ.init = 300;
equ.k = 0;
equ.Q = '4.8*10^16*(T-T2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.001:0.2], ...
                'atol',{'0.000010'}, ...
                'rtol',0.0001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.2    Surface: Temperature [K]', ...
         'axis',[-0.5875993250835259,0.5875993250835259,-0.5155318592631256,0.21553184734219663]);

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5);

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5);

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '1.61*10000/(300*10^(-15))';
equ.init = 300;
equ.k = '386*10^8';
equ.Q = '4.8*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '3.4*10^6/(300*10^(-15))';
equ.init = 300;
equ.k = 0;
equ.Q = '4.8*10^16*(T-T2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.001:0.2], ...
                'atol',{'0.000010'}, ...
                'rtol',0.0001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.2    Surface: Temperature [K]', ...
         'axis',[-0.587599325083526,0.5875993250835262,-0.5155318592631257,0.21553184734219671]);

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '1.61*10000/(300*10^(-15))';
equ.init = 300;
equ.k = '386*10^8';
equ.Q = '4.8*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '3.4*10^6/(300*10^(-15))';
equ.init = 300;
equ.k = 0;
equ.Q = '4.8*10^16*(T-T2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.001:0.2], ...
                'atol',{'0.000010'}, ...
                'rtol',0.0001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.2    Surface: Temperature [K]', ...
         'axis',[-0.587599325083526,0.5875993250835262,-0.5155318592631257,0.21553184734219671]);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '1.61*10000/(300*10^(-15))';
equ.init = 300;
equ.k = '158500*T^(-1.23)*10^8';
equ.Q = '4.8*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '3.4*10^6/(300*10^(-15))';
equ.init = 300;
equ.k = 0;
equ.Q = '4.8*10^16*(T-T2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.001:0.2], ...
                'atol',{'0.000010'}, ...
                'rtol',0.0001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.2    Surface: Temperature [K]', ...
         'axis',[-0.5993152625468479,0.5993152625468481,-0.5155318592631257,0.21553184734219671]);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '158500*T^(-1.23)*10^8';
equ.Q = '4.8*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '3.4*10^6/(300*10^(-15))';
equ.init = 300;
equ.k = 0;
equ.Q = '4.8*10^16*(T-T2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.001:0.2], ...
                'atol',{'0.000010'}, ...
                'rtol',0.0001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.2    Surface: Temperature [K]', ...
         'axis',[-0.5993152625468476,0.5993152625468481,-0.5155318592631257,0.21553184734219671]);

% Animate solution
postmovie(fem, ...
          'tridata',{'T','cont','internal','unit','K'}, ...
          'trimap','jet(1024)', ...
          'axis',[-0.587599325083526,0.5875993250835264,-0.5405349294614421,0.2405349175405131], ...
          'fps',10);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '158500*T^(-1.23)*10^8';
equ.Q = '4.8*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '3.4*10^6/(300*10^(-15))';
equ.init = 300;
equ.k = 0;
equ.Q = '4.8*10^16*(T-T2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.587599325083526,0.5875993250835264,-0.5405349294614421,0.2405349175405131]);

% Animate solution
postmovie(fem, ...
          'tridata',{'T','cont','internal','unit','K'}, ...
          'trimap','jet(1024)', ...
          'axis',[-0.587599325083526,0.5875993250835264,-0.5405349294614421,0.2405349175405131], ...
          'fps',10);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '158500*T^(-1.23)*10^8';
equ.Q = '4.8*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '158500*T^(-1.23)*10^8';
equ.Q = '4.8*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.5993152625468476,0.5993152625468481,-0.5155318592631257,0.21553184734219671]);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '4.8*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '4.8*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.5993152625468479,0.5993152625468483,-0.5155318592631258,0.2155318473421968]);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '4.8*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 23;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '4.8*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.5993152625468479,0.5993152625468483,-0.5155318592631258,0.2155318473421968]);

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.5993152625468479,0.5993152625468483,-0.5155318592631258,0.2155318473421968]);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '96.6*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '4.8*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '4.8*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.5993152625468479,0.5993152625468483,-0.5155318592631258,0.2155318473421968]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.587599325083526,0.5875993250835264,-0.5358154741925482,0.23581546227161912]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.587599325083526,0.5875993250835264,-0.5381608569781535,0.23816084505722454]);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '96.6*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '4.8*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)*100';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '4.8*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.5993152625468479,0.5993152625468483,-0.5155318592631258,0.2155318473421968]);

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.587599325083526,0.5875993250835264,-0.5405349294614422,0.24053491754051318]);

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.587599325083526,0.5875993250835264,-0.5381608569781535,0.23816084505722454]);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '96.6*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '4.8*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '4.8*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.6364171008250145,0.6364171008250149,-0.5381608569781535,0.23816084505722454]);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '50*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '4.8*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '4.8*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.6364171008250145,0.6364171008250149,-0.5381608569781535,0.23816084505722454]);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '70*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '4.8*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '4.8*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.6364171008250145,0.6364171008250149,-0.5381608569781535,0.23816084505722454]);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '70*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '136*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '136*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.6364171008250145,0.6364171008250149,-0.5381608569781535,0.23816084505722454]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.587599325083526,0.5875993250835264,-0.5381608569781535,0.23816084505722454]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.587599325083526,0.5875993250835264,-0.5405349294614421,0.24053491754051318]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'triz',{'T','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]   Height: Temperature [K]', ...
         'grid','on');

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '70*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '3.6*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '3.6*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {'g','3*10^16*(T/645)*(4/(645/T2)*int(x^4/(e^x-1),0,645/T2)-(T2/T)*4/(645/T2)*int(x^4/(e^x-1),0,645/T2))'};

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'triz',{'T','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]   Height: Temperature [K]', ...
         'grid','on');

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.5,0.5,-0.30000001192092896,0]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.587599325083526,0.5875993250835264,-0.5381608569781535,0.23816084505722454]);

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.587599325083526,0.5875993250835264,-0.5381608569781535,0.23816084505722454]);

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.587599325083526,0.5875993250835264,-0.5381608569781535,0.23816084505722454]);

% Plot in cross-section or along domain
postcrossplot(fem,1,[0 1;0 0], ...
              'lindata','T', ...
              'title','Temperature [K]', ...
              'axislabel',{'Arc-length','Temperature [K]'});

% Plot in cross-section or along domain
postcrossplot(fem,1,[0 1;0 0], ...
              'surfdata','T', ...
              'surfmap','jet(1024)', ...
              'title','Temperature [K]');

% Plot in cross-section or along domain
postcrossplot(fem,0,[], ...
              'pointdata',{'T','unit','K'}, ...
              'title','Temperature [K]', ...
              'axislabel',{'Time','Temperature [K]'}, ...
              'refine','auto');

% Plot in cross-section or along domain
postcrossplot(fem,1,[4], ...
              'lindata','T', ...
              'cont','internal', ...
              'title','Temperature [K]', ...
              'axislabel',{'Arc-length','Temperature [K]'}, ...
              'refine','auto');

% Plot in cross-section or along domain
postcrossplot(fem,0,[3], ...
              'pointdata',{'T','unit','K'}, ...
              'title','Temperature [K]', ...
              'axislabel',{'Time','Temperature [K]'}, ...
              'refine','auto');

% Plot in cross-section or along domain
postcrossplot(fem,0,[4], ...
              'pointdata',{'T','unit','K'}, ...
              'pointxdata',{'T','unit','K'}, ...
              'title','Temperature [K] versus Temperature [K]', ...
              'axislabel',{'Temperature [K]','Temperature [K]'}, ...
              'refine','auto');

% Plot in cross-section or along domain
postcrossplot(fem,0,[4], ...
              'pointdata',{'T','unit','K'}, ...
              'pointxdata',{'T','unit','K'}, ...
              'title','Temperature [K] versus Temperature [K]', ...
              'axislabel',{'Temperature [K]','Temperature [K]'}, ...
              'refine','auto');

% Plot in cross-section or along domain
postcrossplot(fem,1,[3], ...
              'lindata','T', ...
              'cont','internal', ...
              'title','Temperature [K]', ...
              'axislabel',{'Arc-length','Temperature [K]'}, ...
              'refine','auto');

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.6364171008250145,0.6364171008250151,-0.5381608569781535,0.23816084505722454]);

% Plot in cross-section or along domain
postcrossplot(fem,1,[3], ...
              'lindata','T', ...
              'cont','internal', ...
              'title','Temperature [K]', ...
              'axislabel',{'Arc-length','Temperature [K]'}, ...
              'refine','auto');

% Plot in cross-section or along domain
postcrossplot(fem,1,[2], ...
              'lindata','T', ...
              'cont','internal', ...
              'title','Temperature [K]', ...
              'axislabel',{'Arc-length','Temperature [K]'}, ...
              'refine','auto');

% Plot in cross-section or along domain
postcrossplot(fem,1,[1], ...
              'lindata','T', ...
              'cont','internal', ...
              'title','Temperature [K]', ...
              'axislabel',{'Arc-length','Temperature [K]'}, ...
              'refine','auto');

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '70*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '3.6*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '3.6*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {'g','3*10^16*(T/645)*(4/(645/T2)*int(x^4/(e^x-1),0,645/T2)-(T2/T)*4/(645/T2)*int(x^4/(e^x-1),0,645/T2))'};

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.6239758642675329,0.6239758642675336,-0.5381608569781535,0.23816084505722454]);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5);

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '70*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '3.6*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '3.6*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {'g','3*10^16*(T/645)*(4/(645/T2)*int(x^4/(e^x-1),0,645/T2)-(T2/T)*4/(645/T2)*int(x^4/(e^x-1),0,645/T2))'};

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.6239758642675329,0.6239758642675336,-0.5381608569781535,0.23816084505722454]);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5);

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '70*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '3.6*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '3.6*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {'g','3*10^16*(T/645)*(4/(645/T2)*int(x^4/(e^x-1),0,645/T2)-(T2/T)*4/(645/T2)*int(x^4/(e^x-1),0,645/T2))'};

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.6239758642675329,0.6239758642675336,-0.5381608569781535,0.23816084505722454]);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5);

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '70*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '3.6*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '3.6*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {'g','3*10^16*(T/645)*(4/(645/T2)*int(x^4/(e^x-1),0,645/T2)-(T2/T)*4/(645/T2)*int(x^4/(e^x-1),0,645/T2))'};

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.6239758642675329,0.6239758642675336,-0.5381608569781535,0.23816084505722454]);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',8);

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '70*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '3.6*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '3.6*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {'g','3*10^16*(T/645)*(4/(645/T2)*int(x^4/(e^x-1),0,645/T2)-(T2/T)*4/(645/T2)*int(x^4/(e^x-1),0,645/T2))'};

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.6626043675817191,0.6626043675817198,-0.5621907491240224,0.2621907372030934]);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',8);

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'boxcoord',[-0.23615740322603818 -0.11095059101695748 -0.10103298534404187 -0.03409072931146401], ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hmax',[100], ...
                  'hmaxfact',1.9, ...
                  'hcurve',0.6, ...
                  'hgrad',1.5, ...
                  'hcutoff',0.01);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '70*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '3.6*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '3.6*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {'g','3*10^16*(T/645)*(4/(645/T2)*int(x^4/(e^x-1),0,645/T2)-(T2/T)*4/(645/T2)*int(x^4/(e^x-1),0,645/T2))'};

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0.001034    Surface: Temperature [K]', ...
         'axis',[-0.6626043675817191,0.6626043675817198,-0.5621907491240224,0.2621907372030934]);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',4);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',4);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',3);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',2);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',3);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',9);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',8);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',7);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '70*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '3.6*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '3.6*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {'g','3*10^16*(T/645)*(4/(645/T2)*int(x^4/(e^x-1),0,645/T2)-(T2/T)*4/(645/T2)*int(x^4/(e^x-1),0,645/T2))'};

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.6626043675817191,0.6626043675817198,-0.5621907491240224,0.2621907372030934]);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',6);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '70*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '3.6*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '3.6*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {'g','3*10^16*(T/645)*(4/(645/T2)*int(x^4/(e^x-1),0,645/T2)-(T2/T)*4/(645/T2)*int(x^4/(e^x-1),0,645/T2))'};

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.6626043675817191,0.6626043675817198,-0.5621907491240224,0.2621907372030934]);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',4);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '70*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '3.6*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '3.6*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {'g','3*10^16*(T/645)*(4/(645/T2)*int(x^4/(e^x-1),0,645/T2)-(T2/T)*4/(645/T2)*int(x^4/(e^x-1),0,645/T2))'};

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.6626043675817191,0.6626043675817198,-0.5621907491240224,0.2621907372030934]);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '70*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '3.6*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '3.6*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {'g','3*10^16*(T/645)*(4/(645/T2)*int(x^4/(e^x-1),0,645/T2)-(T2/T)*4/(645/T2)*int(x^4/(e^x-1),0,645/T2))'};

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.6626043675817191,0.6626043675817198,-0.5621907491240224,0.2621907372030934]);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',3);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '70*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '3.6*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '3.6*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {'g','3*10^16*(T/645)*(4/(645/T2)*int(x^4/(e^x-1),0,645/T2)-(T2/T)*4/(645/T2)*int(x^4/(e^x-1),0,645/T2))'};

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.6626043675817191,0.6626043675817198,-0.5621907491240224,0.2621907372030934]);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',2);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '70*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '3.6*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '3.6*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {'g','3*10^16*(T/645)*(4/(645/T2)*int(x^4/(e^x-1),0,645/T2)-(T2/T)*4/(645/T2)*int(x^4/(e^x-1),0,645/T2))'};

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.6626043675817191,0.6626043675817198,-0.5621907491240224,0.2621907372030934]);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hmaxfact',0.3, ...
                  'hcurve',0.25, ...
                  'hgrad',1.2, ...
                  'hcutoff',0.0003);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',2, ...
                  'point',[], ...
                  'edge',[], ...
                  'subdomain',[], ...
                  'meshstart',fem.mesh);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',2);

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '70*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '3.6*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '3.6*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {'g','3*10^16*(T/645)*(4/(645/T2)*int(x^4/(e^x-1),0,645/T2)-(T2/T)*4/(645/T2)*int(x^4/(e^x-1),0,645/T2))'};

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=0 (2)    Surface: Temperature [K]', ...
         'axis',[-0.6626043675817191,0.6626043675817198,-0.5621907491240224,0.2621907372030934]);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',2);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5);

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '70*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '3.6*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '3.6*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {'g','3*10^16*(T/645)*(4/(645/T2)*int(x^4/(e^x-1),0,645/T2)-(T2/T)*4/(645/T2)*int(x^4/(e^x-1),0,645/T2))'};

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.6626043675817191,0.6626043675817198,-0.5621907491240224,0.2621907372030934]);

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '70*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '3.6*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '3.6*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {'g','3*10^16*(T/645)*(4/(645/T2)*int(x^4/(e^x-1),0,645/T2)-(T2/T)*4/(645/T2)*int(x^4/(e^x-1),0,645/T2))'};

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.6626043675817191,0.6626043675817198,-0.5621907491240224,0.2621907372030934]);

% Plot solution
postplot(fem, ...
         'contdata',{'T','cont','internal','unit','K'}, ...
         'contlevels',20, ...
         'contlabel','off', ...
         'contmap','cool(1024)', ...
         'solnum','end', ...
         'title','Time=2    Contour: Temperature [K]', ...
         'axis',[-0.6239758642675329,0.6239758642675336,-0.5621907491240224,0.2621907372030934]);

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.6239758642675329,0.6239758642675336,-0.5621907491240224,0.2621907372030934]);

% Create mapped quad mesh
fem.mesh=meshmap(fem, ...
                 'hauto',4);

% Create mapped quad mesh
fem.mesh=meshmap(fem, ...
                 'hauto',4);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',4);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',4, ...
                  'jiggle','off');

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',4, ...
                  'jiggle','off');

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',4, ...
                  'jiggle','off');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '70*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '3.6*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '3.6*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {'g','3*10^16*(T/645)*(4/(645/T2)*int(x^4/(e^x-1),0,645/T2)-(T2/T)*4/(645/T2)*int(x^4/(e^x-1),0,645/T2))'};

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.6626043675817191,0.6626043675817198,-0.5621907491240224,0.2621907372030934]);

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',3, ...
                  'jiggle','off');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '70*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '3.6*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '3.6*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {'g','3*10^16*(T/645)*(4/(645/T2)*int(x^4/(e^x-1),0,645/T2)-(T2/T)*4/(645/T2)*int(x^4/(e^x-1),0,645/T2))'};

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.6626043675817191,0.6626043675817198,-0.5621907491240224,0.2621907372030934]);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',2, ...
                  'jiggle','off');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',8, ...
                  'jiggle','off');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',8, ...
                  'jiggle','off');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '70*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '3.6*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '3.6*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {'g','3*10^16*(T/645)*(4/(645/T2)*int(x^4/(e^x-1),0,645/T2)-(T2/T)*4/(645/T2)*int(x^4/(e^x-1),0,645/T2))'};

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=3.689834e-11    Surface: Temperature [K]', ...
         'axis',[-0.6626043675817191,0.6626043675817198,-0.5621907491240224,0.2621907372030934]);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '70*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '3.6*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '3.6*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {'g','3*10^16*(T/645)*(4/(645/T2)*int(x^4/(e^x-1),0,645/T2)-(T2/T)*4/(645/T2)*int(x^4/(e^x-1),0,645/T2))'};

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=6.420267e-12    Surface: Temperature [K]', ...
         'axis',[-0.6239758642675329,0.6239758642675336,-0.5621907491240224,0.2621907372030934]);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',8, ...
                  'jiggle','off');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',7, ...
                  'jiggle','off');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',6, ...
                  'jiggle','off');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',7, ...
                  'jiggle','off');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',8, ...
                  'jiggle','off');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',9, ...
                  'jiggle','off');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '70*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '3.6*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '3.6*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {'g','3*10^16*(T/645)*(4/(645/T2)*int(x^4/(e^x-1),0,645/T2)-(T2/T)*4/(645/T2)*int(x^4/(e^x-1),0,645/T2))'};

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=1.865028e-13    Surface: Temperature [K]', ...
         'axis',[-0.6626043675817191,0.6626043675817198,-0.5621907491240224,0.2621907372030934]);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',5, ...
                  'jiggle','off');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');
% COMSOL Multiphysics Model M-file
% Generated by COMSOL 3.4 (COMSOL 3.4.0.248, $Date: 2007/10/10 16:07:51 $)

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=1.865028e-13    Surface: Temperature [K]', ...
         'axis',[-0.6642919814378613,0.6986877862341371,-0.573577972433116,0.2735779605121869]);

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=1.865028e-13    Surface: Temperature [K]', ...
         'axis',[-0.6642919814378613,0.6986877862341371,-0.573577972433116,0.2735779605121869]);

% Initialize mesh
fem.mesh=meshinit(fem, ...
                  'hauto',4, ...
                  'jiggle','off');

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '70*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '3.6*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '3.6*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {'g','3*10^16*(T/645)*(4/(645/T2)*int(x^4/(e^x-1),0,645/T2)-(T2/T)*4/(645/T2)*int(x^4/(e^x-1),0,645/T2))'};

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.6637116406985003,0.6981074454947761,-0.6005520916910759,0.30055207977014686]);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '80*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '3.6*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '3.6*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {'g','3*10^16*(T/645)*(4/(645/T2)*int(x^4/(e^x-1),0,645/T2)-(T2/T)*4/(645/T2)*int(x^4/(e^x-1),0,645/T2))'};

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2    Surface: Temperature [K]', ...
         'axis',[-0.7202816331340872,0.754677437930363,-0.5998004336190341,0.29980042169810495]);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '80*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '3.6*10^16*(T2-T)+2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '3.6*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {'g','3*10^16*(T/645)*(4/(645/T2)*int(x^4/(e^x-1),0,645/T2)-(T2/T)*4/(645/T2)*int(x^4/(e^x-1),0,645/T2))'};

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2.4], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2.4    Surface: Temperature [K]', ...
         'axis',[-0.7202816331340872,0.754677437930363,-0.5998004336190341,0.29980042169810495]);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '80*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '3.6*10^16*(T2-T)+5*2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '3.6*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {'g','3*10^16*(T/645)*(4/(645/T2)*int(x^4/(e^x-1),0,645/T2)-(T2/T)*4/(645/T2)*int(x^4/(e^x-1),0,645/T2))'};

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2.4], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2.4    Surface: Temperature [K]', ...
         'axis',[-0.7202816331340872,0.754677437930363,-0.5998004336190341,0.29980042169810495]);

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% Refine mesh
fem.mesh=meshrefine(fem, ...
                    'mcase',0, ...
                    'rmethod','regular');

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '80*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '3.6*10^16*(T2-T)+5*2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '3.6*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {'g','3*10^16*(T/645)*(4/(645/T2)*int(x^4/(e^x-1),0,645/T2)-(T2/T)*4/(645/T2)*int(x^4/(e^x-1),0,645/T2))'};

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2.4], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2.4    Surface: Temperature [K]', ...
         'axis',[-0.6637116406985003,0.6981074454947761,-0.600552091691076,0.3005520797701468]);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '80*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '3.6*10^16*(T2-T)+6*2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '3.6*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {'g','3*10^16*(T/645)*(4/(645/T2)*int(x^4/(e^x-1),0,645/T2)-(T2/T)*4/(645/T2)*int(x^4/(e^x-1),0,645/T2))'};

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2.4], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2.4    Surface: Temperature [K]', ...
         'axis',[-0.7202816331340872,0.754677437930363,-0.5998004336190341,0.29980042169810495]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2.4    Surface: Temperature [K]', ...
         'axis',[-0.6637116406985003,0.6981074454947761,-0.5998004336190341,0.29980042169810495]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2.4    Surface: Temperature [K]', ...
         'axis',[-0.6637116406985003,0.6981074454947761,-0.5998004336190341,0.29980042169810495]);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '80*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '1.6*10^16*(T2-T)+6*2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '1.6*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {'g','3*10^16*(T/645)*(4/(645/T2)*int(x^4/(e^x-1),0,645/T2)-(T2/T)*4/(645/T2)*int(x^4/(e^x-1),0,645/T2))'};

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2.4], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2.4    Surface: Temperature [K]', ...
         'axis',[-0.7202816331340872,0.754677437930363,-0.5998004336190341,0.29980042169810495]);

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2.4    Surface: Temperature [K]', ...
         'axis',[-0.6637116406985003,0.6981074454947761,-0.5998004336190341,0.29980042169810495]);

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2.4    Surface: Temperature [K]', ...
         'axis',[-0.6637116406985003,0.6981074454947761,-0.5998004336190341,0.29980042169810495]);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '80*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '0.6*10^16*(T2-T)+6*2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '0.6*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {'g','3*10^16*(T/645)*(4/(645/T2)*int(x^4/(e^x-1),0,645/T2)-(T2/T)*4/(645/T2)*int(x^4/(e^x-1),0,645/T2))'};

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2.4], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2.4    Surface: Temperature [K]', ...
         'axis',[-0.7202816331340872,0.754677437930363,-0.5998004336190341,0.29980042169810495]);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '80*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '0.6*10^16*(T2-T)+7*2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '0.6*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {'g','3*10^16*(T/645)*(4/(645/T2)*int(x^4/(e^x-1),0,645/T2)-(T2/T)*4/(645/T2)*int(x^4/(e^x-1),0,645/T2))'};

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2.4], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2.4    Surface: Temperature [K]', ...
         'axis',[-0.7202816331340872,0.754677437930363,-0.5998004336190341,0.29980042169810495]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2.4    Surface: Temperature [K]', ...
         'axis',[-0.6637116406985003,0.6981074454947761,-0.5998004336190341,0.29980042169810495]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2.4    Surface: Temperature [K]', ...
         'axis',[-0.6637116406985003,0.6981074454947761,-0.5998004336190341,0.29980042169810495]);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '80*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '0.6*10^16*(T2-T)+9*2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2)';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '0.6*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {'g','3*10^16*(T/645)*(4/(645/T2)*int(x^4/(e^x-1),0,645/T2)-(T2/T)*4/(645/T2)*int(x^4/(e^x-1),0,645/T2))'};

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2.4], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2.4    Surface: Temperature [K]', ...
         'axis',[-0.7202816331340872,0.754677437930363,-0.5998004336190341,0.29980042169810495]);

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2.4    Surface: Temperature [K]', ...
         'axis',[-0.6637116406985003,0.6981074454947761,-0.5998004336190341,0.29980042169810495]);

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2.4    Surface: Temperature [K]', ...
         'axis',[-0.6637116406985003,0.6981074454947761,-0.5998004336190341,0.29980042169810495]);

% Animate solution
postmovie(fem, ...
          'tridata',{'T','cont','internal','unit','K'}, ...
          'trimap','jet(1024)', ...
          'axis',[-0.6637116406985003,0.6981074454947761,-0.5998004336190341,0.29980042169810495], ...
          'fps',10);

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2.4    Surface: Temperature [K]', ...
         'axis',[-0.6637116406985003,0.6981074454947761,-0.5998004336190341,0.29980042169810495]);

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum',25, ...
         'title','Time=0.24    Surface: Temperature [K]', ...
         'axis',[-0.6637116406985003,0.6981074454947761,-0.5998004336190341,0.29980042169810495]);

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum',25, ...
         'title','Time=0.24    Surface: Temperature [K]', ...
         'axis',[-0.6637116406985003,0.6981074454947761,-0.5998004336190341,0.29980042169810495]);

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum',62, ...
         'title','Time=0.61    Surface: Temperature [K]', ...
         'axis',[-0.6637116406985003,0.6981074454947761,-0.5998004336190341,0.29980042169810495]);

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum',118, ...
         'title','Time=1.17    Surface: Temperature [K]', ...
         'axis',[-0.6637116406985003,0.6981074454947761,-0.5998004336190341,0.29980042169810495]);

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum',118, ...
         'title','Time=1.17    Surface: Temperature [K]', ...
         'axis',[-0.6637116406985003,0.6981074454947761,-0.5998004336190341,0.29980042169810495]);

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2.4    Surface: Temperature [K]', ...
         'axis',[-0.6637116406985003,0.6981074454947761,-0.5998004336190341,0.29980042169810495]);

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2.4    Surface: Temperature [K]', ...
         'axis',[-0.6637116406985003,0.6981074454947761,-0.5998004336190341,0.29980042169810495]);

% Animate solution
postmovie(fem, ...
          'tridata',{'T','cont','internal','unit','K'}, ...
          'trimap','jet(1024)', ...
          'axis',[-0.6637116406985003,0.6981074454947761,-0.5998004336190341,0.29980042169810495], ...
          'fps',10);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2.4    Surface: Temperature [K]', ...
         'axis',[-0.6637116406985003,0.6981074454947761,-0.5998004336190341,0.29980042169810495]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2.4    Surface: Temperature [K]', ...
         'axis',[-0.6637116406985003,0.6981074454947761,-0.5998004336190341,0.29980042169810495]);

% Animate solution
postmovie(fem, ...
          'tridata',{'T2','cont','internal','unit','K'}, ...
          'trimap','jet(1024)', ...
          'axis',[-0.6637116406985003,0.6981074454947761,-0.5998004336190341,0.29980042169810495], ...
          'fps',10);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum',25, ...
         'title','Time=0.24    Surface: Temperature [K]', ...
         'axis',[-0.6637116406985003,0.6981074454947761,-0.5998004336190341,0.29980042169810495]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum',25, ...
         'title','Time=0.24    Surface: Temperature [K]', ...
         'axis',[-0.6637116406985003,0.6981074454947761,-0.5998004336190341,0.29980042169810495]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum',62, ...
         'title','Time=0.61    Surface: Temperature [K]', ...
         'axis',[-0.6637116406985003,0.6981074454947761,-0.5998004336190341,0.29980042169810495]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum',62, ...
         'title','Time=0.61    Surface: Temperature [K]', ...
         'axis',[-0.6637116406985003,0.6981074454947761,-0.5998004336190341,0.29980042169810495]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum',118, ...
         'title','Time=1.17    Surface: Temperature [K]', ...
         'axis',[-0.6637116406985003,0.6981074454947761,-0.5998004336190341,0.29980042169810495]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum',118, ...
         'title','Time=1.17    Surface: Temperature [K]', ...
         'axis',[-0.6637116406985003,0.6981074454947761,-0.5998004336190341,0.29980042169810495]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2.4    Surface: Temperature [K]', ...
         'axis',[-0.6637116406985003,0.6981074454947761,-0.5998004336190341,0.29980042169810495]);

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2.4    Surface: Temperature [K]', ...
         'axis',[-0.6637116406985003,0.6981074454947761,-0.5998004336190341,0.29980042169810495]);

% (Default values are not included)

% Application mode 1
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.assignsuffix = '_ht';
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '80*T/(300*10^(-15))';
equ.init = 300;
equ.k = '(-0.556+7.13*10^(-3)*T)*10^8';
equ.Q = '0.6*10^16*(T2-T)+9*2*10^23*exp(-1*(-1*z*30+14*t^2))*exp(-50*r^2/(2.5^2))';
equ.rho = 1;
equ.ind = [1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2
clear appl
appl.mode.class = 'HeatTransfer';
appl.mode.type = 'axi';
appl.dim = {'T2'};
appl.name = 'ht2';
appl.assignsuffix = '_ht2';
clear prop
clear weakconstr
weakconstr.value = 'off';
weakconstr.dim = {'lm2'};
prop.weakconstr = weakconstr;
appl.prop = prop;
clear bnd
bnd.type = {'T','q0'};
bnd.T0 = {300,0};
bnd.ind = [1,1,2,1];
appl.bnd = bnd;
clear equ
equ.C = '880/(300*10^(-15))';
equ.init = 300;
equ.k = '(158500*T2^(-1.23))*10^8';
equ.Q = '0.6*10^16*(T-T2)';
equ.rho = 2329;
equ.ind = [1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'r','z'};
fem.frame = {'ref'};
fem.border = 1;
fem.outform = 'general';
clear units;
units.basesystem = 'SI';
fem.units = units;

% Scalar expressions
fem.expr = {'g','3*10^16*(T/645)*(4/(645/T2)*int(x^4/(e^x-1),0,645/T2)-(T2/T)*4/(645/T2)*int(x^4/(e^x-1),0,645/T2))'};

% ODE Settings
clear ode
clear units;
units.basesystem = 'SI';
ode.units = units;
fem.ode=ode;
% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

% Solve problem
fem.sol=femtime(fem, ...
                'solcomp',{'T','T2'}, ...
                'outcomp',{'T','T2'}, ...
                'tlist',[0:0.01:2.4], ...
                'atol',{'0.00010'}, ...
                'rtol',0.001, ...
                'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

% Plot solution
postplot(fem, ...
         'tridata',{'T2','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2.4    Surface: Temperature [K]', ...
         'axis',[-0.7202816331340872,0.7546774379303632,-0.5998004336190341,0.29980042169810495]);

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2.4    Surface: Temperature [K]', ...
         'axis',[-0.7058647399582693,0.7402605447545453,-0.624760169924591,0.32476015800366187]);

% Plot solution
postplot(fem, ...
         'tridata',{'T','cont','internal','unit','K'}, ...
         'trimap','jet(1024)', ...
         'solnum','end', ...
         'title','Time=2.4    Surface: Temperature [K]', ...
         'axis',[-0.7058647399582693,0.7402605447545453,-0.627646249948689,0.32764623802775983]);