I have found different values of F_11 in terms of T_f, but I am not gating values of T_f in eular method.

clear
clc
n = 172;                       %no. of day
H = 20000;                     % altitude in meter 
CF = 0.5;                      %cloud fraction
V_wind =10;                    %wind velocity   
%%Constants taken
u = [7.447323    2.071808    9.010095    7.981491 194];
% a = u(1); % 7.447323;
% b = u(2); % 2.071808;
% c = u(3); % 9.010095;
% d = u(4); % 7.981491;
l = u(5); % 194;

%%constants
alpha = 0.15;                                              %absorptivity of film
albedo_g = 0.35;                                           %ground albedo                                     
albedo_c = 0.5;                                            %cloud albedo
P_sea = 101325;                                            %pressure at sea level
e = 0.0167;                                                %eccentricity of earth
T_g = 280;                                                 %temperature of ground in k
T_cl = 262.15;                                             % assuming cloud at altitude 4km in k
R = 6371000;                                               %redius of earth
epsilon_e = 0.95;                                          %earth emissivity
alpha_ir = 0.85;                                           %infrared absorbtivity of film 
sigma = 5.67.*10.^-8;                                      %stefan boltzman constant 
epsilon_sky = 1;                                                   % emissivity of sky
Pw = 0.0042;
g = 9.81;
[T_oa, P_oa, rho_oa]  = atmos(H);
x = (0:0.1:u(5));
Fun = 1/8*(sqrt(u(1)*(u(5)-x).*(u(2).*x-u(5)*sqrt(u(3))+sqrt(u(3)*u(5).^2-u(4)*u(5).*x))));
D = 2*max(Fun);
syms x 
Fun = 1/8*(sqrt(u(1)*(u(5)-x).*(u(2).*x-u(5)*sqrt(u(3))+sqrt(u(3)*u(5).^2-u(4)*u(5).*x))));
h=1;
k = 10;
zi = k:k:360;
x=0:h:l;                                                             
Fun1 = matlabFunction(Fun);
area = 0;
Q_D = 0;
Q_d=0;
Q_r=0;
Q_earth = 0;
Q_sky = 0;
Q_IR = 0;
Q_aIR = 0;
Q_oc = 0;
% Q_sun = 0;
t = 6:1:18;
Q_sun = zeros(1,13);
for c = 1:(length(t));
   w = 15.*(12-t(c)); % hour angle   
   phi = 18.975;                                                      % lattitude
   delta = 23.45.*sin((pi/180).*(360/365).*(284+n));                  %declination angle
   omega = asin(sin((pi/180).*delta).*sin((pi/180).*phi)+cos((pi/180).*delta).*cos((pi/180).*phi).*cos((pi/180).*w)); %elevation angle
   psi = asin((sin(pi/180).*w).*cos((pi/180).*delta)/cos(omega));     % azimuth angle 
   E = [cos(omega)*cos(psi), -sin(psi)*cos(omega), sin(omega)];       % solar vector
   MA = 2.*pi.*n/365;
   M = (P_oa/P_sea).*(sqrt(1229+(614.*sin(omega)).^2)-614.*sin(omega));
   Tau_atm = 0.5.*(exp(-0.65.*M)+exp(-0.95.*M));
   zita = MA+2.*e.*sin(M)+1.25.*e.^2.*sin(2.*M); 
   ID = (1367.*Tau_atm.^M).*((1+e.*cos(zita))./(1-e.^2)).^2;
   for j=1:length(zi)
       for i=1:(length(x)-1)
          y1 = Fun1(x(i));
          z1 = Fun1(x(i));
          x2 = x(i+1);
          y2 = Fun1(x(i+1));
          z2 = Fun1(x(i+1));
          x3 = x(i+1);
          y3 = y2*sin((pi/180)*zi(j));
          z3 = z2*cos((pi/180)*zi(j)); 
          c1 = y2*cos((pi/180)*k);
          c2 = z2*sin((pi/180)*k);
          T1 = [0, y3-y2, z3-z2];
          T2 = [x(i)-x(i+1), y1-y2, z1-z2];
          vec = cross(T1,T2);                                                % normal vector
          deltal = sqrt((y2-y1)^2+h^2);
          deltazi = sqrt(c2^2+(y2-c1)^2);
          dA11 = deltal.*deltazi;                                            
          area = area +dA11;
          beta = acos((dot(E,vec)/(norm(E).*norm(vec))));                    % angle between two vector
          thetai = (pi/2)+asin(vec(3)/(sqrt(vec(1)^2+vec(2)^2))); 
          QD = alpha.*dA11.*ID.*sin(omega).*cos(pi-beta);                   % direct radiation in Watt
          Q_D = nansum(QD);
          I_d = (1/2).*ID.*sin(omega).*(1.-Tau_atm.^M)./(1.-1.4.*log(Tau_atm));
          Qd = alpha.*I_d.*dA11.*(1+cos(thetai))./2;                         % diffuse radiation
          Q_d1 = nansum(Qd);
          Q_d = Q_d + Q_d1;

            %Reflected solar radiation model
            albedo = albedo_g.*(1-CF).^2+albedo_c.*CF;
            Ir = albedo.*(ID.*sin(omega)+I_d);
            Qr = alpha.*Ir.*dA11.*(1-cos(thetai))/2;                          % reflected radiation
            Q_r = nansum(Qr);
            if beta <= pi/2                                                    % this loop is for consideration of direct radiation,  
                Q_sun(c) = (Q_sun(c)+Q_D+Q_d+Q_r);                                       
            else Q_sun(c) = (Q_sun(c)+Q_d+Q_r);
            end

            %Infrared Radiation from earth
            tau_atm_IR = 1.716-0.5.*(exp(-0.65.*P_oa/P_sea) + exp(-0.95.*P_oa/P_sea));
            T_e = T_g.*(1-CF)+T_cl.*CF;                                        %temp of earth
            eta = asin(R/(R+H));
            VF = (1-cos(eta))/2;                                               %view factor
            Q_e = sigma.*epsilon_e.*alpha_ir.*VF.*(T_e.^4).*tau_atm_IR.*dA11;  %radiation from earth
            Q_earth = Q_earth+Q_e;

            %Infrared Radiation from sky

            T_sky = T_oa.*(0.51+0.208.*sqrt(Pw)).^0.25;                        % temperature of sky 
            Q_s = sigma.*epsilon_sky.*alpha_ir.*(1-VF).*(T_sky.^4).*dA11;      % radiation from sky
            Q_sky = Q_sky + Q_s;

            %Infrared Radiation from airship
            syms T_f;
            epsilon_f = 0.85;                                                  % emissivity of film
            Q_ir = 2.*sigma.*epsilon_f.*(T_f.^4).*dA11;                        % radiation emits from film 
            Q_IR = Q_IR+Q_ir;                                                  

            r = 0.05;
            Q_a = alpha_ir.*sigma.*epsilon_f.*(T_f.^4/(1-r)).*dA11;            % internal radiation absorped by airship

            Q_aIR = Q_aIR + Q_a;
            Q_rr = epsilon_f.*(dA11/2).*sigma.*((T_f^4-T_sky^4).*(1-VF)+(T_f^4-T_e^4).*VF);
            Q_r = Q_r + Q_rr;

            % Convection heat transfer model
            % external convection
            K_air = 0.0241.*(T_oa/273.15).^0.9;
            beta_air = 1/T_oa;
            Pr_air = 0.804-3.25.*10.^(-4).*T_oa;
            nu_air= 1.46.*10.^-6.*T_oa.^1.5/(rho_oa.*(T_oa+110.4));
            Ra_air = g.*beta_air.*(T_f-T_oa).*D.^3/nu_air.^2.*Pr_air;
            Re_air = V_wind.*l/nu_air;                                                                 %reynold no. 
            h_free = (K_air/D).*(0.6+0.387.*((Ra_air)/(1+(0.559/Pr_air).^(9/16)).^(16/9)).^(1/6)).^2;  %natural convective heat transfer
            h_forced = (K_air/D).*(Re_air.*Pr_air.*(0.2275/(log(Re_air)).^2.584-850/Re_air));          %forced convective heat transfer
            h_oc = (h_free.^3 + h_forced.^3).^(1/3);
            Qoc = h_oc.*(T_f-T_oa).*dA11;                                                                   %outside convective heat transfer 
            Q_oc = Qoc+Q_oc;                                                                                
        end
    end
end
Area1 = area;
% disp(Q_sun);
C_p = 1.42*10^3;                                                             %specific heat of material in J
mass = 120*10.^-3*Area1;
% disp(-Q_r-Q_IR-Q_oc+Q_aIR);
% disp(Q_sky+Q_earth-Q_IR-Q_IR-Q_oc+Q_aIR);
% disp(Q_d)

T_f = zeros(1,length(t));

F_11 = (Q_sun+Q_sky+Q_earth+Q_aIR-Q_IR-Q_oc)/(mass.*C_p);
F = matlabFunction(F_11);
T_f(1) = 216.65;
for a = 1:(length(t)-1);
      T_f(a+1) = T_f(a) + F(a);
end