How to Design a Ray Tracing Program for a Concave Mirror

I wrote a ray tracing code for a concave mirror with a planar light source parallel beam. I set up a coordinate system for the light source plane, then transformed the coordinates to the mirror coordinate system, with the mirror coordinate system as the main coordinate system. I set a light source point A, with the light ray making a 30-degree angle with the mirror normal. I assumed the light ray incident at point P on the mirror surface, set the P point coordinates according to the mirror equation, calculated the incident and reflected ray vectors, and based on Fermat's principle (optical path extremum), established a system of equations to calculate the P coordinates. Then, based on the reflection relationship, I calculated the coordinates of point B on the screen. However, my equation system cannot calculate the coordinates of point P. I would like to ask whether there is an error in how I set up the equation system or if the coordinate points and positions I set are unreasonable.

clc
clear
% r=300;
r=300;
r0_prime=350;
alpha=pi/6;
beta0=-alpha;
syms xi w l L M N L_prime N_prime M_prime bianliang
R=200;
P=[xi,w,l];
% 后续建立循环，设n个A坐标
A=[0,10,20];

s = A(2);
z = A(3);
new_x = r * cos(alpha) + s * sin(alpha);
new_y = r * sin(alpha) - s * cos(alpha);
new_z = z; 
A=[new_x,new_y,new_z];
% disp(new_x);
% P坐标关系
xi = R - sqrt(R^2 - w^2 - l^2);
     % 定义偏导数函数（消去 ξ）
dxi_dw = w / sqrt(R^2 - w^2 - l^2);  % ∂ξ/∂w
dxi_dl = l / sqrt(R^2 - w^2 - l^2);  % ∂ξ/∂l  
AP = sqrt((new_x - xi)^2 + (new_y - w)^2 + (new_z - l)^2);
        L = (new_x - xi) / AP;
        M = (new_y - w) / AP;
        N = (new_z - l) / AP;
        
        T1 = (2 * (L - M * (dxi_dw) - N * (dxi_dl))) / (1 + (dxi_dw)^2 + (dxi_dl)^2);
        L_prime = -L + T1;
        M_prime = -M - T1 * dxi_dw;
        N_prime = -N - T1 * dxi_dl;
eq1 = xi == R - sqrt(R^2 - w^2 - l^2);
eq2 = (xi - 260.3076) / sqrt((xi - 260.3076)^2 + (w - 149.1340)^2 + (l - 2)^2) == -sqrt(3) / 2;
eq3 = -(L + L_prime) * dxi_dw - (M + M_prime) == 0;
% 联立方程组
eqns = [eq1; eq2;eq3]

eqns = 

sol = solve(subs(eqns, R, 1000), [w, l], 'Real', true);

Warning: Possibly spurious solutions.

Warning: Possibly spurious solutions.

Warning: Unable to solve symbolically. Returning a numeric solution using vpasolve.

%% 

        T2 = (r0_prime - xi*cos(beta0) - w*sin(beta0)) / (L_prime*cos(beta0) + M_prime*sin(beta0));
        x_prime = xi + T2 * L_prime;
        y_prime = w + T2 * M_prime; 
        z_prime = l + T2 * N_prime;
        Y = (y_prime - r0_prime * sin(beta0)) * sec(beta0);
        Z = z_prime;
        figure;
plot(Y, Z, 'ro', 'MarkerSize', 10, 'MarkerFaceColor', 'r');

Error using plot
Data must be numeric, datetime, duration, categorical, or an array convertible to double.

grid on;

[Torsten]

If you look at "eqns", you only have one "real" equation, but two unknowns (see above).