definite integration, I don't know if my integration is correct

I have this code, I have to integrate the expression expr, and I use polyfit and polyval to aquire the value of G that give to me the value T_Ninf = T. I think that my code is wrong, in particular the integration that I use is wrong because T_Ninf is a vector with high values.

clear all
close all
clc

%%% Temperature
T0 = 15; %°C
L = 0.0065; %°C/m (lapse rate)
h = 4000; %m
T_4000 = T0-(h*L);
T0k = 273.15+15;
T_4000k = T_4000+273.15; %K

%Pressione
P_4000_P0 = (T_4000k/T0k)^5.256;
P0 = 101325; %N/m^2
R = 287.3; %J/(kg K)
P_isa4000 = P_4000_P0*P0*exp((4000-0)/(287*T_4000k)); %N/m^2

%Densità
rho_4000 = P_isa4000/(R*T_4000k);

%Viscosità dinamica e cinematica
a = 1.46*10^-6; %coefficiente
b = 110; %K
mu = a*(T_4000k^1.5/(b+T_4000k));
nu = mu/rho_4000;

%Velocità del suono
gamma = 1.4;
a = sqrt(R*gamma*T_4000k);

%Stima di massima del numero di Reynold
R = (457*10^-3)/2; %m Il diametro 
c = 0.2*R; %m valutazione al 75% del raggio
Omega = 561.82; %rad/s  5365 rpm
Re_075 = (rho_4000*Omega*0.75*R*0.2*R)/mu;

%Parametri
MM = 3.5; %kg; 
M = MM/4;
% M_correct = M(8);
T = 9.81*M;
% T_correct = Tn(8);
A = R^2*pi; %area disco
r = linspace(10^-4,R,200); %valore estratto dal CAD

%Distribuzione radiale di circuitazione costante
% gam =     (T*2*pi)/(A*rho_4000*Omega); %circuitazione
% u =        gam./(2*pi*r); %velocità di rotazione al di sotto del disco
% u_inf =    u/2;
% vi =       sqrt(T/(2*rho_4000*A));

%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%Teoria dei vortici%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%% numero infinito di pale %%%
G =        linspace (10^-5,1,200);
   
r_s =      r./(R.*G);
theta =    acos((r_s.^6+3*r_s.^4+3*r_s.^2-1)./(r_s.^6+4*r_s.^4+3*r_s.^2+1));
omega_s =  6./(5+r_s.^2+2.*(1+r_s.^2).*cos(theta./3));
gam_opt =  2*pi*Omega.*omega_s.*r.^2; %distribusione ottima di circuitazione sul disco attuatore
u_r =      gam_opt./(2*pi*r);
omega =    u_r./r;
gam_opt_ad = omega_s.*(r./(R*G)).^2; %la velocità indotta è costante per un disco uniformemente caricato

% determinazione G di ottimo (ATTENZIONE !!!)
expr=@(r)2*pi*rho_4000.*(Omega-(omega./2)).*omega.*r.^3; %la funzione precedente è uguale, ma con omega esplicitato
T_Ninf =  integral(expr,0,R,'ArrayValued',true);
p = polyfit(T_Ninf,G,2);
G_opt = polyval(p,T);

[Fabio Freschi]

It seems correct: omega is an array. I tried to numerically calculate the integral, adding this piece of code at the bottom of yours

% sampling of r
rr = linspace(0,R,10000);

% loop over omega
for i = 1:length(omega)
    expr1 = @(r)2*pi*rho_4000.*(Omega-(omega(i)./2)).*omega(i).*r.^3;
    TT(i) = trapz(rr,expr1(rr));
end

In the end TT is practically identical to T_Ninf

[Antonino Alesci]

therefore, I think that somethink is wrong from the physical point of view. Anyway, thanks very much