how to solve equations containing symbolic variable

The following error occurred converting from sym to double:
Error using mupadmex
Error in MuPAD command: DOUBLE cannot convert the input expression into a double array.
If the input expression contains a symbolic variable, use the VPA function instead.

I use vpa but I didn't gate answer programme is

clear all
n = 172
H = 20000
CF = 0.5
V_wind =10
%%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;
%%Profile
x = (0:0.1:u(5));
f = 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(f);
R = D/2;
clear x

%%Area of Surface
syms x 
y = 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))));   % Given curve for shape
y1 = diff(y);
y2 = sqrt(1+y1.^2);
y3 = y*y2;
A = 2*pi*int(y3,0,u(5));
vpa(A,5)                                                        % Area of Envelope
A = ans

m = 120*10.^-3*A*0.12*10.^-3

clear x y
%%constants
alpha = 0.15;
albedo_ground = 0.35;
albedo_cloud = 0.5;
P_sea = 101325;
e = 0.0167;
%%solar model
syms t
delta = 23.45.*sin((pi/180).*(360/365).*(284+n));
% [T_oa, P_oa, rho_oa]  = atmos(H);
P_oa = 5446;
rho_oa = 0.088;
T_oa = 216.65;
w = 15.*(12-t);
phi = 18.975;

theta = asin(sin((pi/180).*delta).*sin((pi/180).*phi)+cos((pi/180).*delta).*cos((pi/180).*phi).*cos((pi/180).*w));
% psi = asin((sin(pi/180).*w).*cos((pi/180).*delta)/cos(theta));
MA = 2.*pi.*n/365;
M = (P_oa/P_sea).*(sqrt(1229+(614.*sin(theta)).^2)-614.*sin(theta));  %call fun
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;
I_d = (1/2).*ID.*sin(theta).*(1.-Tau_atm.^M)./(1.-1.4.*log(Tau_atm));
Q_D = alpha.*ID.*sin(theta)                                  % direct radiation

% diffuse soalr radiation model

Q_d = alpha.*I_d

%Reflected solar radiation model
albedo = albedo_ground.*(1-CF).^2+albedo_cloud.*CF;
Ir = albedo.*(ID.*sin(theta)+I_d);
omega = pi.*abs(12-t)/12; 
Q_r = alpha.*Ir.*cos(omega)

%Infrared Radiation from 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;      
alpha_ir = 0.85;
sigma = 5.67.*10.^-8;
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;
eta = asin(R/(R+H));
VF = (1-cos(eta))/2;
Q_earth = sigma.*epsilon_e.*alpha_ir.*VF.*tau_atm_IR.*T_e.^4

%Infrared Radiation from sky
Pw = 0.0042;  
epsilon_sky = 1;
T_sky = T_oa.*(0.51+0.208.*sqrt(Pw)).^0.25;
Q_sky = sigma.*epsilon_sky.*alpha_ir.*(1-VF).*T_sky.^4

%Infrared Radiation from airship
syms T_f
epsilon_f = 0.85;
Q_a = 2.*sigma.*epsilon_f.*T_f.^4

r = 0.05;
Q_aIR = alpha_ir.*sigma.*epsilon_f.*T_f.^4/(1-r)

% Convection heat transfer model
% external convection
g = 9.81;
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;
h_free = (K_air/D).*(0.6+0.387.*((Ra_air)/(1+(0.559/Pr_air).^(9/16)).^(16/9)).^(1/6)).^2;
h_forced = (K_air/l).*(Re_air.*Pr_air.*(0.2275/(log(Re_air)).^2.584-850/Re_air));
h_oc = (h_free.^3 + h_forced.^3).^(1/3);
Q_oc = h_oc.*(T_f-T_oa)

% %internal convection
% syms T_he
% rho_he = 0.1786
% k_he = 0.144.*(T_he/273.15).^0.7
% Pr_he = 0.729-1.6.*10.^(-4).*T_he
% beta_he = 1/T_he
% nu_he = 4.32.*10.^-6.*(T_he.^0.674)/rho_he
% Ra_he = g.*beta_he.*(T_f-T_he).*D.^3/nu_he.^2.*Pr_he
% if Ra_he<1.5.*10.^8
%     Nu_ic = 2.5.*(2+0.6.*Ra_he.^(1/4))
% else Nu_ic = 0.325.*Ra_he.^(1/3)
% end
%     
% h_ic = Nu_ic.*k_he/D 
% Q_ic=h_ic.*(T_f-T_he)

%%solution
clear t T_f
C_p = 1.42;
h = 1;
T_f = zeros(1,length(t));
T_f(1) = 216;
F_xy = @(t,T_f)(Q_D+Q_d+Q_r+Q_earth+Q_sky-Q_aIR-Q_a-Q_oc)/(m.*C_p)
vpa(F_xy,5)
for i=1:(length(t))                              % calculation loop

      k_1 = F_xy(t(i),T_f(i))
      vpa(k_1, 5)
      k_2 = F_xy(t(i)+0.5.*h,T_f(i)+0.5.*h.*k_1)
      vpa(k_2,5)
      k_3 = F_xy((t(i)+0.5.*h),(T_f(i)+0.5.*h.*k_2))
      vpa(k_3,5)
      k_4 = F_xy((t(i)+h),(T_f(i)+k_3.*h))
      vpa(k_4,5)

      T_f(i+1) = T_f(i) + (1/6).*(k_1+2.*k_2+2.*k_3+k_4).*h  % main equation
  end

Answer 1

[Walter Roberson]

F_xy = @(t,T_f)(Q_D+Q_d+Q_r+Q_earth+Q_sky-Q_aIR-Q_a-Q_oc)/(m.*C_p)

defines an anonymous function that takes two parameters, ignores them, and instead returns the value of that expression. The expression involves a symbolic variable, so whenever you call F_xy the result is going to be a symbolic expression. You have assigned numeric values to T_f but you have not told MATLAB to replace the symbolic variable T_f with those numeric values. You pass numeric t into F_xy but because you are ignoring that t, no association is made between that t and the symbolic t.

When you assign a numeric value to a variable that was earlier named in a symbolic expression, uses of that expression continue to get the symbolic variable, not the numeric variable! You should have a look at subs()

In the meantime, though, replace

F_xy = @(t,T_f)(Q_D+Q_d+Q_r+Q_earth+Q_sky-Q_aIR-Q_a-Q_oc)/(m.*C_p);

with

F_xy = matlabFunction((Q_D+Q_d+Q_r+Q_earth+Q_sky-Q_aIR-Q_a-Q_oc)/(m.*C_p), 'Vars', [t T_f]);