Error: Limits of integration must be double or single scalars. While solving equations numerically with symbolic integral limits

Question

Problem Description:

The equation system to be solved {f1=0, f2=0, f3=0} contains variables pl and pr. The expressions for pl and pr involve integrals with variables as their limits (indicated by bold parts in the code). After processing through matlabFunction to create anonymous functions, "int" is replaced with "integral". However, "integral" does not support variables as limits, resulting in an error: Error using integral Limits of integration must be double or single scalars.

Potentially useful reference:

https://www.zhihu.com/question/46915508 This is a similar problem I found, but I can't understand it... https://ww2.mathworks.cn/matlabcentral/answers/2005242-error-limits-of-integration-must-be-double-or-single-scalars-while-solving-equations-numerically-w#answer_1283977

Problem Code:

++++++++++++++++++++++++++

alph=pi*15/180;
L=0.116;
T_ice=18;
addweight=0.212;
T0=15;
dL=0.02;
tilt=0;
heightofweight=12.5;
forcefromcable=0;
mu_L=0.001;
rho_S=920*0.95;
rho_L=1000;
Cp_s=2049.41;
h_m=334000+Cp_s*T_ice;
K_L=0.57;
P_e=102770;
g1=addweight*9.8;
g2=0.38*9.8;
heightofmasscenter=0.07;
heightofcable=0.13;

syms x1 y1 x Uratio

d = pi*L*sin(alph)/2; 
M = g1*(dL+heightofweight*tan(tilt))*cos(tilt)+g2*heightofmasscenter*sin(tilt)-forcefromcable*heightofcable*sin(tilt);
G = g1+g2;

beta = atan(x1/y1);
xL = sqrt(x1^2+y1^2)*cos(pi-alph-beta);
DL = sqrt(x1^2+y1^2)*sin(pi-alph-beta);
xR = sqrt(x1^2+y1^2)*cos(beta-alph);
DR = sqrt(x1^2+y1^2)*sin(beta-alph);

LL = DL^2+(xL-x)^2;
RR = DR^2+(xR-x)^2;
C = Uratio*(int(x*RR^2,x,0,L)-int(x*LL^2,x,0,-L))/(int(LL^1.5,x,0,-L)-int(RR^1.5,x,0,L));
P0 = 12*mu_L*rho_S^4*h_m^3*Uratio^3*(Uratio*int(LL^2*x,x,0,-L)+C*int(LL^1.5,x,0,-L))/(rho_L*T0^3*K_L^3)+P_e;

pl = -12*mu_L*rho_S^4*h_m^3*Uratio^3*(Uratio*int(LL^2*x,0,x)+C*int(LL^1.5, 0,x))/(rho_L*T0^3*K_L^3)+P0;
pr = -12*mu_L*rho_S^4*h_m^3*Uratio^3*(Uratio*int(RR^2*x,0,x)+C*int(RR^1.5, 0,x))/(rho_L*T0^3*K_L^3)+P0;

% Perform definite integrals on pl and pr to obtain the target equation system
f1 = int(pl*x,x,-L,0)+int(pr*x,x,0,L)+M/d; 
f2 = int(pl,x,-L,0)*sin(alph+tilt)+int(pr,x,0,L)*sin(alph-tilt)-G/d-2*L*P_e*sin(alph)*cos(tilt); 
f3 = int(pl,x,-L,0)*cos(alph+tilt)-int(pr,x,0,L)*cos(alph-tilt)+2*L*P_e*sin(alph)*sin(tilt); 

% Proceed with solving
% Initial guess for variables
initial_guess = [0.5, 0.1, 2];
% Define solving function
equations = matlabFunction(f1,f2,f3, 'Vars', {x1, y1, Uratio});
% Use fsolve to solve the equation system
result = fsolve(@(vars) equations(vars(1), vars(2), vars(3)), initial_guess);`