We are asked to find the best approximation in least-squares norm to a
linear fractional transformation:
$$ y = \frac{ax+b}{cx+d} $$
by fitting curve $y_f = c_1 + c_2 x^{c_3}$ in parameters $c_1,c_2,c_3$
over the interval $(0,4)$. Since this is a real interval, I'll assume
that a real approximation is intended, and further, that there's no
singularity of $y$ in this interval or its right endpoint $x=4$.
In the special cases $c=0$ and $d=0$, an exact expression $y = y_f$ is
possible, and also if $bc-ad = 0$ (the function $y$ is constant). This also eliminates the case $y$ is singular at $x=0$.
Otherwise the problem is a
nonlinear least-squares fit
as the model $y_f$ has a nonlinear dependence on the exponent $c_3$.
Because of that, solution by an iterative method is to be expected.
However the nonlinear dependence is only on one parameter, so using
standard methods to handle the other two parameters, the iteration
can be efficiently restricted to a one dimensional search for the
optimal exponent $c_3$.
Note that the family of curves we are modelling with is closed under addition
of constants and multiplication by constants. Subtract a constant $\frac{a}{c}$
from $y$ (adding it to $y_f$) and divide by constant $\frac{bc-ad}{c^2}$,
and we reduce to a case $y = (x+t)^{-1}$ where $t=\frac{d}{c}$. Since
the interval $[0,4]$ does not contain a singularity, either $t \gt 0$
or $t \lt -4$. Accordingly $y$ is concave up and decreasing (positive)
in the first alternative, or concave down and decreasing (negative)
in the second alternative.
Let's write down the objective function $e(c_1,c_2,c_3) = ||y_f - y||_2^2$
and its partials with respect to the three fitting parameters:
$$ e(c_1,c_2,c_3) = \int_0^4 (c_1 + c_2 x^{c_3} - (x+t)^{-1})^2 dx $$
$$ \frac{\partial e}{\partial c_1} = 2 \int_0^4 ((c_1 + c_2 x^{c_3} - (x+t)^{-1}) dx $$
$$ \frac{\partial e}{\partial c_2} = 2 \int_0^4 ((c_1 + c_2 x^{c_3} - (x+t)^{-1}) x^{c_3} dx $$
$$ \frac{\partial e}{\partial c_3} = 2 \int_0^4 ((c_1 + c_2 x^{c_3} - (x+t)^{-1}) c_2 x^{c_3} \log x dx $$
Setting the partial derivatives to zero gives us the necessary "stationary"
condition for an optimal fit. For fixed exponent $c_3$ the first two of
these equations are linear in $c_1,c_2$ and give their values by elimination.
A sensible numerical scheme is then to intially guess an exponent $c_3$
and refine that guess based on standard one dimensional search techniques
such as Golden Section.