How to find an analytical expression ( polynomial fit ) for this ?

                           % program omegaxraio ( principal )

                              tic
                          %T = 0.12;
                          %R=1
                          %R=[0.1:0.1:0.9 1:1:10];
                          R=[1:1:2];
                          T=[1:1:2];
                          m=1;
                          I_value=zeros(1,length(T)*length(R));
                          for i=1:length(T)
                          for k=1:length(R) 

                          g1 = linspace(0.0001,0.001,10);
                          g2 = linspace(0.00101,0.01,10);
                          g3 = linspace(0.0101,0.1,10);
                          g4 = linspace(0.101,1,10);
                          gg4 = linspace(1.01,1.2,20);

                          g5 = linspace(1.01,10,100);
                          g6 = linspace(10.1,40,100);
                          if (T(i) >= 0.1 && T(i) < 1.0)
                              g = [g3 g4 gg4];
                               %g=0.000001;
                          elseif (T(i) >= 1.0 && T(i) < 10)
                              g = [g3 g4 g5 ];
                          elseif (T(i) >= 1.0 && T(i) <= 100)
                              g = [g4 g5 g6];
                              %g=1;
                          end
                          qQd = zeros(1,length(g));
                          t = zeros(1,length(g));
                          for j=1:length(g)
                              [qQd(j), t(j)] = qQdxraio(g(j),T(i),R(k));
                          end
                          t1 =sum(t);

                          %bolinha aqui pra nao rodar omega
                          %plot(g,qQd) 

                          tic
                          fu3 = g.^5.*qQd./(exp(g.^2/T(i)));

                          %I_value = (1/T^3)*trapz(g,fu3)
                          I_value(m) = (1/T(i)^3)*trapz(g,fu3)

                          % fazer grafico omega x T 

                          t_total = toc + t1;
                          disp(['tempo total de integração(s): ' num2str(t_total)]);
                          disp(['valor da integral: ' num2str(I_value)]);
                          plot(g,fu3)
                          m=m+1;
                          end
                          end
                          tic 
                          %semilogx(R,I_value)
                          [T R' I_value']
                          toc

                         %auxiliar program qQdxraio
                            function [qQd, t] = qQdxraio(gg,TT,R)
                                global T g
                                tic
                                T=TT;
                                g = gg;
                                f=@(b) anggbsR(b,R);
                                x=fzero(f,1.1);

                              %f1=@(b) 2*(1 + (sqrt(pi*T)/(2*g)).*sin(ang2(b,R)/2)).*b ;  % parte difusa
                               f2=@(b)  2*(1 - cos(ang2(b,R))).*b ; % parte especular
                               f3=@(b) 2*(1+ (2*(3.14)) .*sin(ang2(b,R)/2)).*b;  % caso novo
                               %f4=@(b) 2*(1 - cos(ang2(b,R))).*b; % somente especular

                                %xx1 = linspace(0,x,150); % parte difuso
                                xx2 = linspace(x,11,150); % parte especular
                                 xx3 = linspace(0,x,150);
                                 %xx4 = linspace(0,11,150);
                                %ff1 = feval(f1,xx1);
                               ff2 = feval(f2,xx2);
                                ff3 = feval(f3,xx3);
                               %ff4 = feval(f4,xx4);

                                qQd = trapz(xx3,ff3) + trapz(xx2,ff2);%  isotropico + especular
                               %qQd = trapz(xx4,ff4); % somente especular
                                %qQd = trapz(xx1,ff1) + trapz(xx2,ff2); %difuso + especular

                                 % bolinha no t , plot(xx1,ff1,xx2,ff2)

                                 t = toc;  

                            disp(['tempo decorrido na última integração (s): ' num2str(t)] );
                            end

                            function xx2 = ang2(b,R)
                               n_iter = length(b);
                               xx2 = zeros(1,n_iter); 
                               for i=1:length(b)
                                    xx2(i)=anggbsR(b(i),R);
                                end
                            end

                            function xx=anggbsR(b,R)
                                global g s bb
                                bb = b;
                                s= 0.3192; 

                               syms r
                                func =  @(r)  1 - (bb./r).^2 - (g^-2)*((2/15)*((s/R))^9 *(1./(r - 1).^9 - 1./(r + 1).^9 - 9./(8*r).*(1./(r - 1).^8 - 1./(r + 1).^8)) - ((s/R))^3 *(1./(r-1).^3 - 1./(r+1).^3 - 3./(2*r).*(1./(r-1).^2 - 1./(r+1).^2)));
                                minrootgbsR = fzero(func,1.1);
                                %minrootgbsR = new_raph(func, 2, 50, 1e-6);
                                k=1;
                                n_iter = length(minrootgbsR); 
                                for i=1:n_iter
                                    vec = zeros(1,n_iter);
                                    if(isreal(minrootgbsR(i))==1)
                                        vec(k)=minrootgbsR(i);
                                        k=k+1;
                                    end
                                end
                                fun =@(r) 1./(r.^2.*sqrt(1 - (b./r).^2 - (g^-2)*(2/15*(s/R)^9 *(1./(r - 1).^9 - 1./(r + 1).^9 - 9./(8* r).*(1./(r - 1).^8 - 1./(r + 1).^8)) - (s/R)^3 *(1./(r -1).^3 - 1./(r + 1).^3 - 3./(2* r).* (1./(r - 1).^2 - 1./(r + 1).^2)))));
                                xx = real(pi - 2*b*integral(fun,max(vec),Inf,  'RelTol',0,'AbsTol',1E-5));
                            end

I want to find a polynomial expression envolving I_value, R and T ... any ideas ?