getvalues := proc() local v,u,w,T1,a,b,alpha,epsilon;
    v := DocumentTools[GetProperty](w_value, value);
    u := DocumentTools[GetProperty](w_unit, value);
    w := parse(v)*Unit(parse(u));
    v := DocumentTools[GetProperty](T1_value, value);
    u := DocumentTools[GetProperty](T1_unit, value);
    T1 := parse(v)*Unit(parse(u));
    v := DocumentTools[GetProperty](a_value, value);
    u := DocumentTools[GetProperty](a_unit, value);
    a := parse(v)*Unit(parse(u));
    v := DocumentTools[GetProperty](b_value, value);
    u := DocumentTools[GetProperty](b_unit, value);
    b := parse(v)*Unit(parse(u));
    v := DocumentTools[GetProperty](alpha_value, value);
    alpha := parse(v);
    v := DocumentTools[GetProperty](epsilon_value, value);
    epsilon := parse(v);
    (w,T1,a,b,alpha,epsilon)
end:

plugin:=proc() local v,r,T,m,P,Topt;
   v := getvalues();
   r := eval(Meq,[w=v[1],T1=v[2],a=v[3],b=v[4],alpha=v[5],epsilon=v[6],sigma=ScientificConstants[Constant](Stefan_Boltzmann_constant, units), T2=T2*Unit(K)]);
   r := simplify(rhs(r));
   r := evalf(r);
   r := convert(r,unit_free);
   T := Optimization[NLPSolve](r, T2=0..1000);
   m := T[1];
   DocumentTools[SetProperty](mass_value, value, m);
   Topt := eval(T2, T[2]);
   DocumentTools[SetProperty](temp_value, value, Topt);
   P:=plot(r, T2=(Topt-100)..(Topt+100),gridlines);
   DocumentTools[SetProperty](opt, value, P);
end:

plugin():

Warning, limiting number of major iterations has been reached

Warning, unable to evaluate the function to numeric values in the region; see the plotting command's help page to ensure the calling sequence is correct