How to take input values from users symbolically while running a code?

o2 = [0, 0, 0]; % Origin for ain
ain = [26, 0, 0]; % Initial vector for ain
input_axis = [0, 1, 0]; % Axis of rotation for ain (y-axis)
theta1 = deg2rad(10); % Angle of rotation for ain in radians

% Rotation matrix function
rot_matrix = @(axis, theta) cos(theta) * eye(3) + …
sin(theta) * [0, -axis(3), axis(2); axis(3), 0, -axis(1); -axis(2), axis(1), 0] + …
(1 – cos(theta)) * (axis’ * axis);

% Compute the rotated vector for ain
a_rotated = rot_matrix(input_axis, theta1) * (ain’ – o2′) + o2′;
a_final = a_rotated’;
disp(norm(a_final));

cin=[122.95, -20, 0];
c_rotated= rot_matrix(input_axis, theta1) * (cin’ – o2′) + o2′;
naxis=[sin(theta1),0,cos(theta1)];
syms phi;
c_final_rotated=rot_matrix(naxis,phi)*(c_rotated-a_rotated)+a_rotated;

bin = [29.5, 30, 0];
o4 = [13.5, 30, 0]; % Origin for bin
output_axis = [0, 1, 0]; % Axis of rotation for bin (y-axis)
theta2 = deg2rad(10); % Angle of rotation for bin in radians

% Compute the rotated vector for bin
b1_rotated = rot_matrix(output_axis, theta2) * (bin’ – o4′) + o4′;
b1_final = b1_rotated’;
disp(norm(b1_final-o4));

coupler = c_final_rotated’-b1_final;
coupler = subs(coupler, conj(phi), phi);

%%in this code,i want to take theta1 and theta2 values from the user before from doing parametric substitution(t) given below and want to show eq in terms of theta1 and theta2 and phi
eq=norm(coupler)-106;
disp(eq);
%%after doing this ,put values of theta 1 and theta2 in
coupler = subs(coupler, conj(phi), phi);
%%and continue further as below

syms t;
cos_phi = (1 – t^2) / (1 + t^2);
sin_phi = 2 * t / (1 + t^2);

% Substitute parametric forms into coupler components
coupler_parametric = subs(coupler, [cos(phi), sin(phi)], [cos_phi, sin_phi]);

% Display the parametric coupler
disp(‘Parametric form of coupler:’);
disp(coupler_parametric);

syms targetvalue % it might be 3.5 …
normsq = expand(sum(coupler_parametric.^2) – targetvalue^2);

normpoly = simplify(normsq*(t^2+1)^2);
vpa(expand(normpoly),4);
tsolve = solve(normpoly,t,’maxdegree’,4,’returnconditions’,true);
h=vpa(subs(tsolve.t,targetvalue, 106));
%disp(h);
real_solutions = h(imag(h) == 0);
disp(‘Real roots:’);
disp(real_solutions);

angles_rad = 2 * atan(real_solutions);
angles_deg = rad2deg(angles_rad);

% Display angles in degrees
disp(‘Angles in degrees before adjustment:’);
disp(angles_deg);

phi=double(angles_rad(1));

c1_position = double(rot_matrix(naxis,phi) * (c_rotated – a_rotated) + a_rotated);
p=(c1_position’-a_final)’;
%q=(c1_position’-b1_final)’;

angle=acosd(p(2)/norm(p));
disp(angle);
%%at last i want to display value of angle for both theta1 and theta2 equal to 10 degrees

%%please help someone ,i am getting no real solution if i am taking theta1 and theta2 from user as 10 degrees but no error if i am explicitly defining like in this codeo2 = [0, 0, 0]; % Origin for ain
ain = [26, 0, 0]; % Initial vector for ain
input_axis = [0, 1, 0]; % Axis of rotation for ain (y-axis)
theta1 = deg2rad(10); % Angle of rotation for ain in radians