Imagine the quadratic curve with equation

y=y(x)=ax^2+bx+c

On the concave side of this curve there is a circle of radius R. The circle is as close the the extremum of the quadratic as possible without resulting in the curves crossing each other. Write a function which takes as inputs a,b,c, and R and returns the coordinates of the center of the circle.

For example, if

a=1; b=0; c=10; R=pi; 

then the function returns

T = circ_puzz(a,b,c,R)
T =
     0   20.1196044010894

This can be visualized as follows:

 P = @(x) a*x.^2 + b*x + c;  % Quadratic
 C = @(x) real(-sqrt(R^2-(x-T(1)).^2) + T(2));  % Lower half circle
 x = linspace(-R,R,10000);  % Range of plotted data.
 plot(x,C(x),'r',x,-C(x)+2*T(2),'r',x,P(x),T(1),T(2),'*k')
 ylim([0,30])
 axis equal