Speed up the code

Hello

This code is part of a total code. It calculates the geometrical features of the main problem geometry. However, it takes so much time due to mupadmex. How can I improve it? Similar code is written for the solution of x2 and x3 which I have not copied (due to lack of space).

BTW, if i omit the if imaginary condition, would there be any problem? I mean is it possible that double command generate imaginary numbers?

any clever way of speeding up the code is highly appriciated.

%Slopes of three lines
m1=(P(2,2)-P(1,2))/(P(2,1)-P(1,1)); % line 1
m2=(P(3,2)-P(2,2))/(P(3,1)-P(2,1)); % line 2
m3=(P(4,2)-P(3,2))/(P(4,1)-P(3,1)); % line 3


syms x1    %introducing the variables x1,x2,x3
syms x2
syms x3

int=0; % counting the points of intersection - it has to be 2 at the end, in case there is intersection


Points(2,2)=0;   %for preallocation

% intersection of circle with line 3
 y3=P(3,2)+m3*(x3-P(3,1));
c3 = (x3-xc)^2+(y3-yc)^2-R^2;
x_sol3 = vpasolve(c3);
y_sol3 = subs(y3,'x3',x_sol3); %convert sym to numeric value in y3
x_sol3 = double(x_sol3);
y_sol3 = double(y_sol3);
if imag(x_sol3(1))==0 && imag(x_sol3(2))==0 % This is to guarantee there is real intersection (not imaginary)
    for v=1:2
        if x_sol3(v)>P(3,1) && x_sol3(v)<P(4,1)
            if v==2 && abs(x_sol3(1)-x_sol3(2))<1e-6 % means if two roots are the same do not count it
                int;
            else
                int=int+1;
                Points(int,:)=[x_sol3(v) y_sol3(v)];

            end
        end
    end
end 

[Pooneh Shah Malekpoor]

I have now replaced symbolic solver with fsolve to speed up the code; however, it only gives out one of the solution and not the other one...the is much more quick now, but does not printout one of the solutions. what should i do?

m1=(P(2,2)-P(1,2))/(P(2,1)-P(1,1)); % line 1

m2=(P(3,2)-P(2,2))/(P(3,1)-P(2,1)); % line 2

m3=(P(4,2)-P(3,2))/(P(4,1)-P(3,1)); % line 3

int=0; % counting the points of intersection - it has to be 2 at the end, in case there is intersection

Points(2,2)=0;

% intersection of circle with line 3

F = @(x)[(x(1)-xc)^2+((P(3,2)+m3*(x(1)-P(3,1)))-yc)^2-R^2;(x(2)-xc)^2+((P(3,2)+m3*(x(2)-P(3,1)))-yc)^2-R^2];

opts = optimoptions('fsolve', 'Display', 'off');

x = fsolve(@(x) F(x), zeros(1,2), opts);

y3=P(3,2)+m3*(x-P(3,1));