How to use fprintf

How can I use fprintf to print out

NSO2;
    
    CSO2_bulkslurry;
    
    CHSO3;
    
    CSO3;
    
    CSO2_inter;
    
    CHSO3_inter;
   
     CSO3_inter

on the script below:

global S_total C_total Ca_total CSO2_bulkslurry  CHSO3 CSO3 CCO2_bulkslurry CHCO3 CCO3 CCaCO3 CCaSO3 Trial CH P CCO2_in nloop Delt t CH_trial ...
CSO2_inter  CHSO3_inter  CSO3_inter CCO2_inter  CHCO3_inter  CCO3_inter pHtrial pi pHInter V_Headspace F CSO2_in R T HSO2 kga kLa DCa2 DSO2 kLa_CO2 HCO2 DCO2 KSPCaSO3  BETCaSO3 ktot BETCaCO3 MWCaCO3 Kad ...
rhoCaCO3 KCO2 KHCO3 KSO2 KHSO3 Kw MWCaSO3 rhoCaSO3 kCaS03 PercentageAbsorption DHSO3 DSO3 Enhancement NSO2

format shorte

%% Parameters

%Symbol     Value           Description                             Units
R           = 8.314;        % Universal gas constant                (J/molK)
HSO2 		= 149;          % SO2 Henry’s law constant              (m3Pa/mol)
HCO2		= 5.15e+3;      % CO2 Henry’s law constant              (m3Pa/mol)
DCa2 		= 1.39e-9;      % Diffusivity of Ca2+ ions              (m2/s)
DSO2 		= 2.89e-9;      % Diffusivity of SO2                    (m2/s)
DCO2 		= 3.53e-9;      % Diffusivity of CO2                    (m2/s)
KCO2		= 1.7e-3;       % CO2 dissociation equilibrium constant	(mol/m3)
KHCO3		= 6.55e-8;      % HCO3- dissociation equilibrium constant(mol/m3)
KSO2		= 6.24;         % SO2 dissociation equilibrium constant	(mol/m3)
KHSO3		= 5.68e-5;      % HSO3- dissociation equilibrium constant(mol/m3)
Kw          = 5.3E-8 ;      % water dissociation constant            [mol/m3.]^2
KSPCaSO3	= 3.1e-4;       % Solubility product of CaSO3            (mol2/m6)
BETCaSO3	= 10;           % BET specific surface area of CaSO3     (m2/g)
ktot		= 8.825e-3;     % Total dissolution rate constant        (m/s)
BETCaCO3	= 12.54;        % BET specific surface area of CaSO3    (m2/g)
MWCaCO3     = 100.0869;     % Molecular weight of CaCO3             (g/mol)
Kad         = 8.4e-3;       % Adsorption constant                   (m3/mol)
rhoCaCO3	= 2703;         % Density of limestone                  (kg/m3)
MWCaSO3     = 258.30;       % Molecular weight of CaSO3.0.5H2O      (g/mol)
rhoCaSO3	= 2540;         % Density of CaSO3.0.5H2O               (kg/m3)
kCaS03      = 2e-08;        % rate constant for CaSO3               (mol/m^3 s)
P           = 1.01E+05;     % operating pressure                    (Pa)
T           = 323.15;       % Reactor temperature                   (K)
V_Headspace = 4.e-4;        % Volume of the headspace               (m3)% changed t0 -4
F           = 1.66667e-5;   % Flue gas inlet flow rate              (m3/s)
y_so2       = 2000* 1E-6;   % 2000 ppm
%CSO2_in     = y_so2 * P/R/T;
CSO2_in     = 7.28E-02;     % Inlet flue gas SO2 concentration       (mol/m3)
y_co2       = 0.0;
CCO2_in     = y_co2 * P/R/T;% CO2  concentration in gas phase,      (molm^3)
kga 		= 4.14e-6;      % Product of gas-side mass coefficient and interfacial surface area	(mol/m3/Pa/s)
kga         = 4.74e-5;      % suggested value
kLa         = 6e-4;         % Product of liquid-side mass coefficient and interfacial surface area for SO2 (1/s)
kLa_CO2     = 9.598e-5;     % Product of gas-side mass coefficient and interfacial surface area for SO2		(1/s)
DHSO3 		= 2.04E-09;		% Diffusivity of HSO3                    (m2/s)
DSO3 		= 1.70E-09;	% Diffusivity of SO3                    (m2/s)

%% Initialization
pHtrial     = 6.8;
Ca_total    = 0 * 4.879E-02;
C_total     = 0 * 4.879E-02;
S_total     = 0.0; 
CCaCO3      = 0.0;  
CCaSO3      = 0.0;

% SO2_g is a fraction (fr) of CSO2_in 
fr          = 0.01;
SO2_g       =  fr * CSO2_in;

%
CO2_g       = CCO2_in;
CO2_g       = 0.2;
CCO2_inter  = CO2_g ;
%
pg          = 200;
CSO2_inter  = pg/HSO2;

% Initial pH (bulk and interface)
CH_trial    = 1.e-7;
pH1         = 6.8;
pH          = fzero (@HionpH, pH1 );
CH          = 10^(3-pH) ;

%% C_Bulk

%  Relationship between CSO2_bulkslurry, S_total,CH+,KSO2  and KHSO3-
    CSO2_bulkslurry = (S_total*CH^2)/(CH^2 + KSO2*CH + KSO2*KHSO3);

% Relationship between HSO3-, S_total,CH+,KSO2  and KHSO3- 
    CHSO3           = (S_total*KSO2*CH)/(CH^2 + KSO2*CH + KSO2*KHSO3);

% Relationship between SO32-, S_total,CH+,KSO2  and KHSO3- 
    CSO3            = (S_total*KSO2*KHSO3)/(CH^2 + KSO2*CH + KSO2*KHSO3);

% Relationship between CCO2_bulkslurry, C_total,CH+,KCO2  and KHCO3- 
    CCO2_bulkslurry = (C_total*CH^2)/(CH^2 + KCO2*CH + KCO2*KHCO3);

% Relationship between HCO3-, C_total,CH+,KCO2  and KHCO3- 
    CHCO3           = (C_total*KCO2*CH)/(CH^2 + KCO2*CH + KCO2*KHCO3);

% Relationship between CO32-, C_total,CH+,KCO2  and KHCO3- 
    CCO3            = (C_total*KCO2*KHCO3)/(CH^2 + KCO2*CH + KCO2*KHCO3);
%%

pH2                 = 6.0;

pHInter             = fzero (@HionStar, pH2); 

CH                  = 10^(3-pH) ;

% set initial conditions
y0                  = [ SO2_g CO2_g  S_total C_total  Ca_total CCaCO3 CCaSO3 ];
y0                  = y0'; % put as a column vector

% set  trial values for Fsolve eo be used in ODES
Trial(1)            = (SO2_g * R *T ) / HSO2;  % assumed so2 at interface
Trial(2)            = 1;                       % SO2_total at interface
Trial(3)            = CH ;                     % hydrogen at interfafe
%% time stepping parametes
Delt                = 1.0   ;                  % time step changed from 0.01;
t0                  = 0;                       % time for further simulation
t                   = t0;
tmax                = 930.90909090909090909090;% maximum time for simulation for each loop.
nloop               = tmax/Delt ;              % number of steps  to be taken.
nstep               = 33 ; 

%  template for writing the solution.
Results             = [ t  y0'    pH  pHInter];
%% time stepping starts.

for k = 1:nstep
   
    % loop for nloop time steps.
    y = SO2_OdeDriver( y0);
    
    % write the results and repeat.
    pH = 3 - log10 (CH) ;
    tf = t;
    
    % intermediate results 
    Results = [Results; tf y' pH pHInter];   
    
    PercentageAbsorption = ((CSO2_in - Results(:,2))/CSO2_in).*100;
    
    Enhancement = 1 + DHSO3.* (CHSO3_inter - CHSO3)/DSO2.* (CSO2_inter - Results(:,2)) + DSO3.* (CSO3_inter - CSO3)/ DSO2.* (CSO2_inter - Results(:,2)); 
    
    %NSO2 = (SO2_g* R* T - HSO2* Results(:,2))./((1/kga) + HSO2/(Enhancement* kLa));
    
    
    y0 =  y;
end

Results; 
Solution = [ Results(:,1) Results(:,2) Results(:,4) Results(:,9)];

PercentageAbsorption = ((CSO2_in - Results(:,2))/CSO2_in).*100

Enhancement

Exp_Time=[0;960;1920;2880;3840;4800;5760;6720;7680;8640;9600;10560;11520;12480;13440;14400;15360;16320;17280;18240;19200;20160;21120;22080;23040;24000;24960;25920;26880;27840;28800;29760;30720];

Exp_SO2_g =[7.28E-02;3.66E-02;4.13E-02;4.44E-02;4.78E-02;5.02E-02;5.30E-02;5.48E-02;5.70E-02;5.88E-02;6.03E-02;6.18E-02;6.35E-02;6.43E-02;6.52E-02;6.62E-02;6.71E-02;6.81E-02;6.85E-02;6.92E-02;6.98E-02;7.02E-02;7.06E-02;7.11E-02;7.12E-02;7.15E-02;7.18E-02;7.20E-02;7.22E-02;7.23E-02;7.25E-02;7.26E-02;7.27E-02];

Exp_S_total = [0;0.35;0.69;0.99;1.26;1.51;1.72;1.91;2.08;2.23;2.37;2.49;2.59;2.68;2.76;2.83;2.89;2.94;2.99;3.03;3.06;3.09;3.11;3.13;3.15;3.16;3.17;3.18;3.19;3.19;3.20;3.20;3.20];

Exp_pH = [6.87;6.65;6.43;5.68;3.45;3;2.85;2.78;2.63;2.61;2.57;2.56;2.55;2.58;2.59;2.55;2.57;2.57;2.52;2.51;2.51;2.53;2.53;2.51;2.5;2.51;2.49;2.49;2.49;2.47;2.47;2.46;2.46];

figure (1)
plot(Exp_Time,Exp_SO2_g,'kd',Results(:,1), Results(:,2),'k-')
legend('Exp-SO_{2,g}','Mod-SO_{2,g}')
ylabel('Conc (mol/m^{3}','FontSize',12,'FontWeight','bold')
xlabel ('Time (sec)','FontSize',12,'FontWeight','bold')
set(gca,'YColor','k');

figure(2)         
plot(Exp_Time,Exp_S_total,'bo', Results(:,1),Results(:,4),'b-')
legend('Exp-S_{total}','Mod-S_{total}','location','best')
ylabel('Conc (mol/m^{3}','FontSize',12,'FontWeight','bold')
xlabel ('Time (sec)','FontSize',12,'FontWeight','bold')
set(gca,'YColor','k');


figure(3)         
plot(Exp_Time,Exp_pH,'rv',Results(:,1),Results(:,9),'r-')
legend('Exp-pH','Mod-pH','location','best')
ylabel('pH ','FontSize',12,'FontWeight','bold')
xlabel ('Time (sec)','FontSize',12,'FontWeight','bold')
set(gca,'YColor','k');


%% pHmodel
function ph =  HionpH (pH)
global V_Headspace F CSO2_in R T HSO2 kga kLa DCa2 DSO2 kLa_CO2 HCO2 DCO2 KSPCaSO3  BETCaSO3 ktot BETCaCO3 MWCaCO3 Kad rhoCaCO3 KCO2 KHCO3 ...
    KSO2 KHSO3 Kw MWCaSO3 rhoCaSO3 S_total C_total Ca_total

% Ph suppied. calculate CH and set up discrepacncy function.

CH =  10^(3-pH);

% Hydrogen balance from electro-neutrality.
ph =  CH + 2  * Ca_total - ((S_total*KSO2*CH)/(CH^2 + KSO2*CH + KSO2*KHSO3))...
 - 2*((S_total*KSO2*KHSO3)/(CH^2 + KSO2*CH + KSO2*KHSO3))...
 -((C_total*KCO2*CH)/(CH^2 + KCO2*CH + KCO2*KHCO3))...
 - 2*((C_total*KCO2*KHCO3)/(CH^2 + KCO2*CH + KCO2*KHCO3))- Kw/CH ;
end


function ph1 =  HionStar (pH)

global V_Headspace F CSO2_in R T HSO2 kga kLa DCa2 DSO2 kLa_CO2 HCO2 DCO2 KSPCaSO3  BETCaSO3 ktot BETCaCO3 MWCaCO3 Kad rhoCaCO3 KCO2 KHCO3 ...
 KSO2 KHSO3 Kw MWCaSO3 rhoCaSO3 S_total C_total Ca_total CSO2_bulkslurry  CHSO3 CSO3 CSO2_inter  CHSO3_inter  CSO3_inter CCO2_inter  CHCO3_inter ...
 CCO3_inter CCO2_bulkslurry CHCO3 CCO3 pg

CH =  10^(3-pH);

% species concentrations.
CSO2_inter ;
CHSO3_inter = KSO2*CSO2_inter  / CH ;
CSO3_inter = KHSO3  *CHSO3_inter  /CH;
Stotal = CHSO3 + CSO3;
Ctotal = CHCO3 + CCO3 ;
CHSO3_inter = Stotal * (   1+  KHSO3  / CH) ;
CSO3_inter = KHSO3  * CHSO3_inter  /CH;
CHCO3_inter = Ctotal * ( 1+  KHCO3  / CH) ;
CCO3_inter = KHCO3  *CHCO3_inter  /CH ;
Ca   = Ca_total;

ph1 = CH + 2*Ca -(CHSO3_inter + 2*CSO3_inter + CHCO3_inter + 2* CCO3_inter +  Kw/CH );

end

function y = SO2_OdeDriver (y0)

global V_Headspace F CSO2_in R T HSO2 kga kLa DCa2 DSO2 kLa_CO2 HCO2 DCO2 KSPCaSO3  BETCaSO3 ktot BETCaCO3 MWCaCO3 Kad rhoCaCO3 KCO2 KHCO3 KSO2 ...
    KHSO3 Kw MWCaSO3 rhoCaSO3 P CCO2_in S_total C_total Ca_total CSO2_bulkslurry  CHSO3 CSO3 CCO2_bulkslurry CHCO3 CCO3 CCaCO3 CCaSO3 CH  ... 
  nloop Delt t CH_trial Tag pHInter

% take nloop steps  starting at y0.
for k = 1:nloop

  pH_trial = 8.0;
  pH = fzero (@HionpH, pH_trial );
  CH = 10^(3-pH);

% report as pH for ease of reference
  pH =3 -log10 (CH) ;

  % calculate the specieswise bulk concentrations by calling
  
  %% C_Bulk

%  Relationship between CSO2_bulkslurry, S_total,CH+,KSO2  and KHSO3-
    CSO2_bulkslurry = (S_total*CH^2)/(CH^2 + KSO2*CH + KSO2*KHSO3);

% Relationship between HSO3-, S_total,CH+,KSO2  and KHSO3- 
    CHSO3 = (S_total*KSO2*CH)/(CH^2 + KSO2*CH + KSO2*KHSO3);

% Relationship between SO32-, S_total,CH+,KSO2  and KHSO3- 
    CSO3 = (S_total*KSO2*KHSO3)/(CH^2 + KSO2*CH + KSO2*KHSO3);

% Relationship between CCO2_bulkslurry, C_total,CH+,KCO2  and KHCO3- 
    CCO2_bulkslurry = (C_total*CH^2)/(CH^2 + KCO2*CH + KCO2*KHCO3);

% Relationship between HCO3-, C_total,CH+,KCO2  and KHCO3- 
    CHCO3 = (C_total*KCO2*CH)/(CH^2 + KCO2*CH + KCO2*KHCO3);

% Relationship between CO32-, C_total,CH+,KCO2  and KHCO3- 
    CCO3 = (C_total*KCO2*KHCO3)/(CH^2 + KCO2*CH + KCO2*KHCO3);
%%
 
% use  these to find dydt expressions calling the odes.m file.
dydt = odes (t, y0);

% integrate with explicit Euler
y = y0 + dydt' * Delt;

% new results are
SO2_g  = y(1);
CO2_g = y(2);
S_total = y(3) ;
C_total = y(4) ;
Ca_total = y(5);
CCaCO3 = y(6);
CCaSO3 = y(7);

% reassign these for the next time step.
y0 = y ;
t = t +Delt;
CH_trial = CH; 
end
end

function dcdt = odes (t, c)
global V_Headspace F CSO2_in R T HSO2 kga kLa DCa2 DSO2 kLa_CO2 HCO2 DCO2 KSPCaSO3  BETCaSO3 ktot BETCaCO3 MWCaCO3 Kad rhoCaCO3 KCO2 KHCO3 KSO2 ...
    KHSO3 Kw MWCaSO3 rhoCaSO3 CSO2_bulkslurry  CHSO3 CSO3 CCO2_bulkslurry CHCO3 CCO3 CCaCO3 CCaSO3 CH CCO2_in S_total C_total Ca_total Trial ...
    CSO2_inter  CHSO3_inter  CSO3_inter CCO2_inter  CHCO3_inter  CCO3_inter pHtrial pi pHInter pg

%% variable definitions for ease of reference
% c(1) = SO2 in the exit gas.
% c(2) = CO2 in the exit gas.
% c(3) =  total sulfur in bulk liquid
% c(4) =  total carbon in bulk liquid
% c(5) =  total calcium in bulk liquid
% c(6) =  solid CaCO3
% c(7) =  solid CaSO3.

% SO2 in gas phase.
S_total  = c(3) ;
C_total = c(4) ;
Ca_total = c(5) ;
% trial value
pH1 = 7; 
pH = fzero (@HionpH, pH1 );
CH = 10^(3-pH) ;


%% C_Bulk

%  Relationship between CSO2_bulkslurry, S_total,CH+,KSO2  and KHSO3-
    CSO2_bulkslurry = (S_total*CH^2)/(CH^2 + KSO2*CH + KSO2*KHSO3);

% Relationship between HSO3-, S_total,CH+,KSO2  and KHSO3- 
    CHSO3 = (S_total*KSO2*CH)/(CH^2 + KSO2*CH + KSO2*KHSO3);

% Relationship between SO32-, S_total,CH+,KSO2  and KHSO3- 
    CSO3 = (S_total*KSO2*KHSO3)/(CH^2 + KSO2*CH + KSO2*KHSO3);

% Relationship between CCO2_bulkslurry, C_total,CH+,KCO2  and KHCO3- 
    CCO2_bulkslurry = (C_total*CH^2)/(CH^2 + KCO2*CH + KCO2*KHCO3);

% Relationship between HCO3-, C_total,CH+,KCO2  and KHCO3- 
    CHCO3 = (C_total*KCO2*CH)/(CH^2 + KCO2*CH + KCO2*KHCO3);

% Relationship between CO32-, C_total,CH+,KCO2  and KHCO3- 
    CCO3 = (C_total*KCO2*KHCO3)/(CH^2 + KCO2*CH + KCO2*KHCO3);
%%

CH;
pH;


% SO2 in gas phase
% bulk gas at exit conditions
p_exit = c(1)*R*T ;
pg = p_exit; 
Trial ;
Tnew = fsolve ( @NSub , Trial);
% Results are
pHInter = 3-log10(Tnew(3) );
pstar = Tnew(1) * HSO2;
NSO2 = kga *(pg-pstar);

% assign new trial values
Trial = Tnew;
dcdt(1) = (1/V_Headspace) * (F * CSO2_in - F * c(1) ) - NSO2  ;

% CO2 in gas phase.
E_CO2 = 1;
NCO2 = 0;
dcdt(2) = (1/V_Headspace) * (F * CCO2_in - F * c(2)) - NCO2 ;



% total sulfur balance
% CaSO3 cystallization sub-model
n = 3;
RCaSO3 =  0 * 1.62E-08 * exp(-5152/T) *( (c(5)*CSO3/ KSPCaSO3) -1)^n ;
% Limestone dissolution submodel
RCaCO3 = ktot*BETCaCO3*MWCaCO3*CCaCO3*CH *(1 - (Kad*CH)/(1+Kad*CH));
% RCaCO3 = NCO2 + RCaCO3;
dcdt(3) = NSO2 -  RCaSO3;
dcdt(4) = NCO2+ RCaCO3;



% total calcium balance in liquud
dcdt(5) =  RCaCO3- RCaSO3;
if ( c(6) < 1.0e-3 )
    c(6) = 0.0; % cut off caco3 here.
    dcdt(6) = 0.0;
    Time_c = t ; % time when caco3 becomes zero.
    Tag = 1;
else
    dcdt(6) = - RCaCO3 ;
end



% CaSO3 balance
dcdt(7) = RCaSO3;
end

       

function discrep = NSub (Tr)
global V_Headspace F CSO2_in R T HSO2 kga kLa DCa2 DSO2 kLa_CO2 HCO2 DCO2 KSPCaSO3  BETCaSO3 ktot BETCaCO3 MWCaCO3 Kad rhoCaCO3 KCO2 KHCO3 KSO2 ...
    KHSO3 Kw MWCaSO3 rhoCaSO3 S_total C_total Ca_total CSO2_bulkslurry  CHSO3 CSO3 CSO2_inter  CHSO3_inter  CSO3_inter CCO2_inter  CHCO3_inter ...
    CCO3_inter CCO2_bulkslurry CHCO3 CCO3 pg pHInter

Tr;
CSO2_inter          = Tr(1) ;
pstar               = CSO2_inter * HSO2 ;
S1Total             = Tr(2) ;                               % CHSO3 + CSO3 at interface
CH                  = Tr(3) ;                                % h-ion  concentration at interface

S_total             = CHSO3 + CSO3 + CSO2_bulkslurry ;  % in bulk for rate-Liq calculation.
CHSO3_inter         = KSO2*CSO2_inter  / CH ;
CSO3_inter          = KHSO3  *CHSO3_inter  /CH;
Ctotal              = CHCO3 + CCO3 ;
CHCO3_inter         = Ctotal * ( 1+  KHCO3  / CH) ;
CCO3_inter          = KHCO3  *CHCO3_inter  /CH ;
S_total_inter       = CSO2_inter + CHSO3_inter +  CSO3_inter;
df                  = S_total_inter- S_total ; % driving force due to reactions
kga ;
pstar;
pg;
Rate_gas = kga* (pg-pstar) ;
Rate_liq = kLa * df; 

% discrepancies to be resolved.
discrep(1) = S1Total - (CHSO3_inter + CSO3_inter) ;
discrep(2)  = Rate_gas - Rate_liq ;

%  minus of H conce  calculated
mCH_cal =  2*Ca_total -(CHSO3_inter + 2*CSO3_inter + CHCO3_inter + 2* CCO3_inter +  Kw/CH ) ;
discrep(3) = CH + mCH_cal ;

end

[Rik]

What do you mean? You want to send the function to you computer printer? That is not what fprintf does. With such a huge number of global variables, I'm not going to spend time understanding your code to see where you might want to print a result to the command window.

[Dursman Mchabe]

I need to print the values per every time step, for:

Line 310

Line 313

Line 316

Line 341

Line 391

Line 397

Line 398

[Rik]

You put your code in your question as code, instead of attaching it as a file, so while it is directly visible, there are no line numbers. Also, I have no indication that you read the documentation for fprintf. What is your actual question? What have you tried?

fprintf('some variable is now %.2f\n',some_variable)

Would that do the trick?

[Dursman Mchabe]

I have tried:

    fprintf('%4.2f/n',NSO2);
    
    fprintf('%4.2f/n',CSO2_bulkslurry);
    
    fprintf('%4.2f/n',CHSO3);
    
    fprintf('%4.2f/n',CSO3);
    
    
    fprintf('%4.2f/n',CSO2_inter);
    
    fprintf('%4.2f/n',CHSO3_inter);
    
    fprintf('%4.2f/n',CSO3_inter);

But it only gives 6 values per species.

[Dursman Mchabe]

I also tried

for m = 1:length(tf)
    fprintf('%4.2f/n',NSO2);
    
    fprintf('%4.2f/n',CSO2_bulkslurry(m));
    
    fprintf('%4.2f/n',CHSO3(m));
    
    fprintf('%4.2f/n',CSO3(m));
    
    
    fprintf('%4.2f/n',CSO2_inter(m));
    
    fprintf('%4.2f/n',CHSO3_inter(m));
    
    fprintf('%4.2f/n',CSO3_inter(m));
end

But it gives error.

[Walter Roberson]

Newline is \n not /n

[Dursman Mchabe]

Ooops, I try

for m = 1:length(tf)
    fprintf('%4.2f\n',NSO2);
    
    fprintf('%4.2f\n',CSO2_bulkslurry(m));
    
    fprintf('%4.2f\n',CHSO3(m));
    
    fprintf('%4.2f\n',CSO3(m));
    
    
    fprintf('%4.2f\n',CSO2_inter(m));
    
    fprintf('%4.2f\n',CHSO3_inter(m));
    
    fprintf('%4.2f\n',CSO3_inter(m));
    end

[Walter Roberson]

What error are you getting?

[Walter Roberson]

I recommend you add the 'Display', 'none' option combination to your fsolve() call.

[Dursman Mchabe]

After reftifying "\n", the error is gone, however, I only get one value per species.

[Dursman Mchabe]

I have used :

options = optimoptions('fsolve','Display','none');
Tnew = fsolve ( @NSub , Trial,options);

But it did not make a difference.