The purpose of this program is to solve 'x' by Newton iteration.
The author knows that there is a problem with 'vpasolve', but there is still a problem after changing to 'solve'.
This program is composed of three equations, which looks more complicated, but it is not difficult to understand because of the large amount of calculation. I have been troubled by various bugs in this program for a long time, and I hope to get your help.
function F=funfqtotceshi7(x)
global R q_tot fq_tot y i A B C TM K E S
q_rad=(T1^4-x^4)*g/(1/c+1/c-1);
k=0.017+(7E-6)*(800-Tm)+0.0228*log(Tm);
syms xn1 xn2 xn3 xn4 xn5 xn6
eq1=vpasolve((C1P*a+C2f*(0.017+(7E-6)*(800-(x+xn1)/2)+0.0228*log((x+xn1)/2))/D)*(x-xn1)+(x.^4-xn1.^4)*g/(1/c+1/c-1)-q_tot==0,xn1);
eq2=vpasolve((C1P*a+C2f*(0.017+(7E-6)*(800-(eq1+xn2)/2)+0.0228*log((eq1+xn2)/2))/D)*(eq1-xn2)+(eq1.^4-xn2.^4)*g/(1/c+1/c-1)-q_tot==0,xn2);
eq3=vpasolve((C1P*a+C2f*(0.017+(7E-6)*(800-(eq2+xn3)/2)+0.0228*log((eq2+xn3)/2))/D)*(eq2-xn3)+(eq2.^4-xn3.^4)*g/(1/c+1/c-1)-q_tot==0,xn3);
eq4=vpasolve((C1P*a+C2f*(0.017+(7E-6)*(800-(eq3+xn4)/2)+0.0228*log((eq3+xn4)/2))/D)*(eq3-xn4)+(eq3.^4-xn4.^4)*g/(1/c+1/c-1)-q_tot==0,xn4);
eq5=vpasolve((C1P*a+C2f*(0.017+(7E-6)*(800-(eq4+xn5)/2)+0.0228*log((eq4+xn5)/2))/D)*(eq4-xn5)+(eq4.^4-xn5.^4)*g/(1/c+1/c-1)-q_tot==0,xn5);
eq6=vpasolve((C1P*a+C2f*(0.017+(7E-6)*(800-(eq5+xn6)/2)+0.0228*log((eq5+xn6)/2))/D)*(eq5-xn6)+(eq5.^4-xn6.^4)*g/(1/c+1/c-1)-q_tot==0,xn6);
A=[T1 x eq1 eq2 eq3 eq4 eq5 eq6];
B(i)=(A(i).^4-A(i+1).^4)*g/(1/c+1/c-1);
C(i)=C1P*a*(A(i)-A(i+1));
K(i)=0.017+(7E-6)*(800-TM(i))+0.0228*log(TM(i));
E(i)=(C2f*K(i)/D)*(A(i)-A(i+1));
R(i)=1/((B(i)+C(i)+E(i))/(A(i)-A(i+1)));
y(i+1)=TC+((sum(R(8-i:7)))/S)*(TH-TC);
fq_tot=(y(8)^4-y(7)^4)*g/(1/c+1/c-1)+C1P*a*(y(8)-y(7))+(C2f*(0.017+(7E-6)*(800-(y(8)+y(7))/2)+0.0228*log((y(8)+y(7))/2))/D)*(y(8)-y(7));
Differentiate Equation dF
function dF=dfunfqtotceshi4(x)
global R q_tot fq_tot y i A B C TM K E S
q_rad=(T1^4-x^4)*g/(1/c+1/c-1);
k=0.017+(7E-6)*(800-Tm)+0.0228*log(Tm);
syms xn1 xn2 xn3 xn4 xn5 xn6
eq1=vpasolve((C1P*a+C2f*(0.017+(7E-6)*(800-(x+xn1)/2)+0.0228*log((x+xn1)/2))/D)*(x-xn1)+(x.^4-xn1.^4)*g/(1/c+1/c-1)-q_tot==0,xn1);
eq2=vpasolve((C1P*a+C2f*(0.017+(7E-6)*(800-(eq1+xn2)/2)+0.0228*log((eq1+xn2)/2))/D)*(eq1-xn2)+(eq1.^4-xn2.^4)*g/(1/c+1/c-1)-q_tot==0,xn2);
eq3=vpasolve((C1P*a+C2f*(0.017+(7E-6)*(800-(eq2+xn3)/2)+0.0228*log((eq2+xn3)/2))/D)*(eq2-xn3)+(eq2.^4-xn3.^4)*g/(1/c+1/c-1)-q_tot==0,xn3);
eq4=vpasolve((C1P*a+C2f*(0.017+(7E-6)*(800-(eq3+xn4)/2)+0.0228*log((eq3+xn4)/2))/D)*(eq3-xn4)+(eq3.^4-xn4.^4)*g/(1/c+1/c-1)-q_tot==0,xn4);
eq5=vpasolve((C1P*a+C2f*(0.017+(7E-6)*(800-(eq4+xn5)/2)+0.0228*log((eq4+xn5)/2))/D)*(eq4-xn5)+(eq4.^4-xn5.^4)*g/(1/c+1/c-1)-q_tot==0,xn5);
eq6=vpasolve((C1P*a+C2f*(0.017+(7E-6)*(800-(eq5+xn6)/2)+0.0228*log((eq5+xn6)/2))/D)*(eq5-xn6)+(eq5.^4-xn6.^4)*g/(1/c+1/c-1)-q_tot==0,xn6);
A=[T1 x eq1 eq2 eq3 eq4 eq5 eq6];
B(i)=(A(i).^4-A(i+1).^4)*g/(1/c+1/c-1);
C(i)=C1P*a*(A(i)-A(i+1));
K(i)=0.017+(7E-6)*(800-TM(i))+0.0228*log(TM(i));
E(i)=(C2f*K(i)/D)*(A(i)-A(i+1));
R(i)=1/((B(i)+C(i)+E(i))/(A(i)-A(i+1)));
y(i+1)=TC+((sum(R(8-i:7)))/S)*(TH-TC);
fq_tot=(y(8)^4-y(7)^4)*g/(1/c+1/c-1)+C1P*a*(y(8)-y(7))+(C2f*(0.017+(7E-6)*(800-(y(8)+y(7))/2)+0.0228*log((y(8)+y(7))/2))/D)*(y(8)-y(7));
Newton iteration
dt(ii)=funfqtotceshi7(x(ii))./dfunfqtotceshi4(x(ii));
disp('number of iterations:');
Many thanks in adavance,
Hopefully there is a solution for this problem.