clear;

model = mphload('heat_transient_axi');

model.param.set('Tinput','1000[degC]');
model.physics('ht').feature('temp1').set('T0', 1, 'Tinput');

pdom1 = model.probe.create('pdom1', 'DomainPoint');
pdom1.model('mod1');
pdom1.setIndex('coords2','0.28',0,0);
pdom1.setIndex('coords2','0.38',0,1);


M = mphstate(model,'sol1','out',{'Mc' 'MA' 'MB' 'C' 'D'},...
'input','T0', 'output', 'mod1.ppb1');

norm(M.MA)
norm(M.MB)
norm(M.Mc)
norm(M.C)
norm(M.D)


T0 = 273.15;
Tinput = 1273.15-T0;

opt = odeset('mass', M.Mc);
func = @(t,x) M.MA*x + M.MB*Tinput;
[t,x] = ode23s(func, [0:10:190], zeros(size(M.MA,1),1), opt);
y = M.C*x';
y = y+T0;

%func1 = @(t,x) M.MA*x;
%[t,xx] = ode23s(func1, [0:10:190], zeros(size(M.MA,1),1), opt);
%yy = M.C*xx';
%yy = yy+T0;

plot(t,y)
hold on
%plot(t,yy,'m*')
%hold on
Tnum = mphinterp(model,'T','coord',[0.28;0.38],'t',t);
plot(t,Tnum,'r+')