clc;clear all;close all;
model = mphopen('COMSOL_BTE_11_19');
info = mphsolinfo(model);
%info_1 = mphsolutioninfo(model, 'dataset', dset1);
% mphmodel(model) Get model nodes
% set a handle for your solution:

solu = mpheval(model,'u');
solu2 = mpheval(model,'u2');
solu3 = mpheval(model,'u3');
solu4 = mpheval(model,'u4');
% which gives a structure that has the solution and the mesh (not the extended mesh!)
soluvec_u = solu.d1;
soluvec_u2 = solu2.d1;
soluvec_u3 = solu3.d1;
soluvec_u4 = solu4.d1;
I_w = [soluvec_u;soluvec_u2;soluvec_u3;soluvec_u4];
n_x = size(I_w);
[Iw_0] = function_gauss_quadrature(I_w);
% I_guess = Iw_0;
c = num2cell(Iw_0);
model.param.set('I_guess',{c{1}}, 'W/(m^2)');
% model.physics('I').feature('cfeq1').set('f', 'I_guess');
% model.physics('I_2').feature('cfeq1').set('f', 'I_guess');
% model.physics('I_3').feature('cfeq1').set('f', 'I_guess');
% model.physics('I_4').feature('cfeq1').set('f', 'I_guess');
model.sol('sol1').runAll;
mphplot(model,'pg1');
new_solu = mpheval(model,'u');
new_solu2 = mpheval(model,'u2');
new_solu3 = mpheval(model,'u3');
new_solu4 = mpheval(model,'u4');