title: Convert OCamCalib toolbox Matlab script to Python
body: <p>I am trying to convert a script from the OCamCalib Matlab toolbox to python. This scripts estimate a inverse polynomial function used for the reverse path of light ray on a fisheye lens.</p>
  <p>My python scripts fails to compute the polynomial values at the end.</p>
  <p>The Matlab script:
  findinvpoly.m</p>
  <pre class="lang-matlab prettyprint-override"><code>%FINDINVPOLY finds the inverse polynomial specified in the argument.
  %   [POL, ERR, N] = FINDINVPOLY(SS, RADIUS, N) finds an approximation of the inverse polynomial specified in OCAM_MODEL.SS.
  %   The returned polynomial POL is used in WORLD2CAM_FAST to compute the reprojected point very efficiently.
  %   
  %   SS is the polynomial which describe the mirrror/lens model.
  %   RADIUS is the radius (pixels) of the omnidirectional picture.
  %   ERR is the error (pixel) that you commit in using the returned
  %   polynomial instead of the inverse SS. N is searched so that
  %   that ERR is &lt; 0.01 pixels.
  %
  %   Copyright (C) 2008 DAVIDE SCARAMUZZA, ETH Zurich
  %   Author: Davide Scaramuzza - email: davide.scaramuzza@ieee.org
  
  function [pol, err, N] = findinvpoly(ss, radius)
  
  if nargin &lt; 3
      maxerr = inf;
      N = 1;
      while maxerr &gt; 0.01 %Repeat until the reprojection error is smaller than 0.01 pixels
          N = N + 1;
          [pol, err] = findinvpoly2(ss, radius, N);
          maxerr = max(err);  
      end
  else
      [pol, err, N] = findinvpoly2(ss, radius, N)
  end
  
  function [pol, err, N] = findinvpoly2(ss, radius, N)
  
  theta = [-pi/2:0.01:1.20];
  r     = invFUN(ss, theta, radius);
  ind   = find(r~=inf);
  theta = theta(ind);
  r     = r(ind);
  
  pol = polyfit(theta,r,N);
  err = abs( r - polyval(pol, theta)); %approximation error in pixels
  
  function r=invFUN(ss, theta, radius)
  
  m=tan(theta);
  
  r=[];
  poly_coef=ss(end:-1:1);
  poly_coef_tmp=poly_coef;
  for j=1:length(m)
      poly_coef_tmp(end-1)=poly_coef(end-1)-m(j);
      rhoTmp=roots(poly_coef_tmp);
      res=rhoTmp(find(imag(rhoTmp)==0 &amp; rhoTmp&gt;0 &amp; rhoTmp&lt;radius ));
      if isempty(res) | length(res)&gt;1
          r(j)=inf;
      else
          r(j)=res;
      end
  end
  </code></pre>
  <p>The script calling the matlab function :
  main.m</p>
  <pre><code>coeff = [-1.729819e+03, 0.000000e+00, 2.218967e-04, -2.655544e-08, 2.169698e-11 ];
  
  
  h = 2160;
  w = 4096;
  radius = sqrt(h*h + w*w);
  
  [pol, err, N] = findinvpoly(coeff, radius);
  
  pol
  
  </code></pre>
  <p>The output I have is :</p>
  <pre><code>pol =
  
     1.0e+03 *
  
      0.0002    0.0020    0.0057    0.0002   -0.0163   -0.0015    0.0345    0.0068   -0.0082    0.0463    0.0192    0.0751    0.2371   -0.0350    1.2974    2.4535
  
  </code></pre>
  <p>My Python script :</p>
  <pre class="lang-py prettyprint-override"><code>import numpy as np
  import math
  
  def find_inverse_polynome(poly_coefficients, radius, N=1, threshold_error = 0.01):
  
      max_error = np.inf
      while max_error &gt; threshold_error:
          N = N+1
          poly, error = compute_inverse_polynome(poly_coefficients, radius, N)
          max_error = error.max()
  
      return poly, error, N
  
  def compute_inverse_polynome(poly_coefficients, radius, N):
      theta = np.arange(-np.pi, 1.2, 0.01)
      r = np.array(inverse_poly(poly_coefficients, theta, radius))
      indices = np.where(r != np.inf)[0]
      theta = theta[indices]
      r = r[indices]
  
      poly = np.polyfit(theta, r, N)
      print(theta)
      exit()
      error = abs(r - np.polyval(poly, theta)) # approximation of error in pixels
  
      return poly, error
  
  def inverse_poly(poly_coefficients, theta, radius):
  
      m = np.tan(theta)
      result = []
      reversed_poly_coefficients = poly_coefficients[::-1]
      reversed_poly_coefficients_temp = reversed_poly_coefficients
  
      for idx, j in enumerate(range(len(m))):
          reversed_poly_coefficients_temp[-2] = reversed_poly_coefficients[-2] - m[j]
          rho_temp = np.roots(reversed_poly_coefficients_temp)
          
          res = rho_temp[np.where((np.imag(rho_temp) == 0) &amp; (rho_temp &lt; radius) &amp; (rho_temp&gt;0))]
          if res.size == 0 or res.size &gt; 1:
              result.append(np.inf)
          else:
              result.append(res)
  
  
      return result
  
  poly_coeff = [-1.729819e+03, 0.000000e+00, 2.218967e-04, -2.655544e-08, 2.169698e-11 ]
  
  poly_coeff = [1707.2, 0,    -0.00023066,    2.4174E-08, -1.8933E-11]
  
  h = 2160
  w = 4096
  radius = math.sqrt(h*h + w*w)
  a = find_inverse_polynome(poly_coeff, radius)
  
  
  </code></pre>
  <p>The current error I am having is :</p>
  <pre><code>Traceback (most recent call last):
    File &quot;/media/alix/System/Users/aleroy/Documents/DPhil/algorithms/opticalflow/generate_invpoly.py&quot;, line 52, in &lt;module&gt;
      a = find_inverse_polynome(poly_coeff, radius)
    File &quot;/media/alix/System/Users/aleroy/Documents/DPhil/algorithms/opticalflow/generate_invpoly.py&quot;, line 9, in find_inverse_polynome
      poly, error = compute_inverse_polynome(poly_coefficients, radius, N)
    File &quot;/media/alix/System/Users/aleroy/Documents/DPhil/algorithms/opticalflow/generate_invpoly.py&quot;, line 22, in compute_inverse_polynome
      poly = np.polyfit(theta, r, N)
    File &quot;&lt;__array_function__ internals&gt;&quot;, line 5, in polyfit
    File &quot;/home/alix/Documents/dphil/venv/p38/lib/python3.8/site-packages/numpy/lib/polynomial.py&quot;, line 629, in polyfit
      c, resids, rank, s = lstsq(lhs, rhs, rcond)
    File &quot;&lt;__array_function__ internals&gt;&quot;, line 5, in lstsq
    File &quot;/home/alix/Documents/dphil/venv/p38/lib/python3.8/site-packages/numpy/linalg/linalg.py&quot;, line 2306, in lstsq
      x, resids, rank, s = gufunc(a, b, rcond, signature=signature, extobj=extobj)
  TypeError: No loop matching the specified signature and casting was found for ufunc lstsq_n
  </code></pre>
  <p><strong>EDIT :</strong>
  I was able to modify the code and it now outputs values.  However the system does  not converge. One condition was missing in the np.where of the inverse_poly() function. Now, I have this warning message :</p>
  <pre><code>RankWarning: Polyfit may be poorly conditioned
    poly, error = compute_inverse_polynome(poly_coefficients, radius, N)
  </code></pre>
  <p>As requested in the comments, this is the current value of r:</p>
  <pre><code>[[2.44779367e+03+0.j]
   [2.43423568e+03+0.j]
   [2.40712477e+03+0.j]
   [2.36649200e+03+0.j]
   [2.31243826e+03+0.j]
   [2.24520354e+03+0.j]
   [2.16526289e+03+0.j]
   [2.07344374e+03+0.j]
   [1.97104699e+03+0.j]
   [1.85993558e+03+0.j]
   [1.74253978e+03+0.j]
   [1.62173327e+03+0.j]
   [1.50057400e+03+0.j]
   [1.38196726e+03+0.j]
   [1.26835405e+03+0.j]
   [1.16151748e+03+0.j]
   [1.06254023e+03+0.j]
   [9.71883778e+02+0.j]
   [8.89530856e+02+0.j]
   [8.15138353e+02+0.j]
   [7.48169221e+02+0.j]
   [6.87991946e+02+0.j]
   [6.33948462e+02+0.j]
   [5.85396629e+02+0.j]
   [5.41734273e+02+0.j]
   [5.02410876e+02+0.j]
   [4.66931443e+02+0.j]
   [4.34855756e+02+0.j]
   [4.05795080e+02+0.j]
   [3.79407661e+02+0.j]
   [3.55393816e+02+0.j]
   [3.33491061e+02+0.j]
   [3.13469549e+02+0.j]
   [2.95127906e+02+0.j]
   [2.78289525e+02+0.j]
   [2.62799306e+02+0.j]
   [2.48520811e+02+0.j]
   [2.35333796e+02+0.j]
   [2.23132080e+02+0.j]
   [2.11821707e+02+0.j]
   [2.01319365e+02+0.j]
   [1.91551022e+02+0.j]
   [1.82450751e+02+0.j]
   [1.73959723e+02+0.j]
   [1.66025327e+02+0.j]
   [1.58600421e+02+0.j]
   [1.51642675e+02+0.j]
   [1.45114005e+02+0.j]
   [1.38980084e+02+0.j]
   [1.33209908e+02+0.j]
   [1.27775427e+02+0.j]
   [1.22651214e+02+0.j]
   [1.17814185e+02+0.j]
   [1.13243344e+02+0.j]
   [1.08919561e+02+0.j]
   [1.04825384e+02+0.j]
   [1.00944862e+02+0.j]
   [9.72633960e+01+0.j]
   [9.37676047e+01+0.j]
   [9.04452052e+01+0.j]
   [8.72849075e+01+0.j]
   [8.42763201e+01+0.j]
   [8.14098661e+01+0.j]
   [7.86767082e+01+0.j]
   [7.60686813e+01+0.j]
   [7.35782325e+01+0.j]
   [7.11983671e+01+0.j]
   [6.89225998e+01+0.j]
   [6.67449111e+01+0.j]
   [6.46597079e+01+0.j]
   [6.26617875e+01+0.j]
   [6.07463058e+01+0.j]
   [5.89087481e+01+0.j]
   [5.71449021e+01+0.j]
   [5.54508347e+01+0.j]
   [5.38228693e+01+0.j]
   [5.22575668e+01+0.j]
   [5.07517067e+01+0.j]
   [4.93022710e+01+0.j]
   [4.79064291e+01+0.j]
   [4.65615237e+01+0.j]
   [4.52650585e+01+0.j]
   [4.40146865e+01+0.j]
   [4.28081993e+01+0.j]
   [4.16435173e+01+0.j]
   [4.05186810e+01+0.j]
   [3.94318424e+01+0.j]
   [3.83812575e+01+0.j]
   [3.73652794e+01+0.j]
   [3.63823519e+01+0.j]
   [3.54310030e+01+0.j]
   [3.45098401e+01+0.j]
   [3.36175442e+01+0.j]
   [3.27528658e+01+0.j]
   [3.19146198e+01+0.j]
   [3.11016821e+01+0.j]
   [3.03129852e+01+0.j]
   [2.95475152e+01+0.j]
   [2.88043081e+01+0.j]
   [2.80824470e+01+0.j]
   [2.73810591e+01+0.j]
   [2.66993132e+01+0.j]
   [2.60364170e+01+0.j]
   [2.53916149e+01+0.j]
   [2.47641859e+01+0.j]
   [2.41534412e+01+0.j]
   [2.35587228e+01+0.j]
   [2.29794012e+01+0.j]
   [2.24148742e+01+0.j]
   [2.18645646e+01+0.j]
   [2.13279196e+01+0.j]
   [2.08044085e+01+0.j]
   [2.02935217e+01+0.j]
   [1.97947697e+01+0.j]
   [1.93076813e+01+0.j]
   [1.88318027e+01+0.j]
   [1.83666965e+01+0.j]
   [1.79119404e+01+0.j]
   [1.74671260e+01+0.j]
   [1.70318581e+01+0.j]
   [1.66057535e+01+0.j]
   [1.61884399e+01+0.j]
   [1.57795551e+01+0.j]
   [1.53787458e+01+0.j]
   [1.49856666e+01+0.j]
   [1.45999792e+01+0.j]
   [1.42213508e+01+0.j]
   [1.38494533e+01+0.j]
   [1.34839622e+01+0.j]
   [1.31245549e+01+0.j]
   [1.27709093e+01+0.j]
   [1.24227025e+01+0.j]
   [1.20796086e+01+0.j]
   [1.17412968e+01+0.j]
   [1.14074293e+01+0.j]
   [1.10776577e+01+0.j]
   [1.07516204e+01+0.j]
   [1.04289382e+01+0.j]
   [1.01092096e+01+0.j]
   [9.79200415e+00+0.j]
   [9.47685546e+00+0.j]
   [9.16325079e+00+0.j]
   [8.85061852e+00+0.j]
   [8.53831107e+00+0.j]
   [8.22558216e+00+0.j]
   [7.91155518e+00+0.j]
   [7.59517858e+00+0.j]
   [7.27516091e+00+0.j]
   [6.94987288e+00+0.j]
   [6.61719378e+00+0.j]
   [6.27425842e+00+0.j]
   [5.91701393e+00+0.j]
   [5.53937971e+00+0.j]
   [5.13147807e+00+0.j]
   [4.67530431e+00+0.j]
   [4.13134483e+00+0.j]
   [3.37490430e+00+0.j]
   [9.69109951e-01+0.j]
   [1.03280934e+00+0.j]
   [1.06639997e+00+0.j]
   [1.08970138e+00+0.j]
   [1.10772769e+00+0.j]
   [1.12251856e+00+0.j]
   [1.13511044e+00+0.j]
   [1.14610442e+00+0.j]
   [1.15588125e+00+0.j]
   [1.16469793e+00+0.j]
   [1.17273640e+00+0.j]
   [1.18013037e+00+0.j]
   [1.18698108e+00+0.j]
   [1.19336713e+00+0.j]
   [1.19935080e+00+0.j]
   [1.20498229e+00+0.j]
   [1.21030271e+00+0.j]
   [1.21534615e+00+0.j]
   [1.22014122e+00+0.j]
   [1.22471214e+00+0.j]
   [1.22907965e+00+0.j]
   [1.23326161e+00+0.j]
   [1.23727355e+00+0.j]
   [1.24112901e+00+0.j]
   [1.24483994e+00+0.j]
   [1.24841688e+00+0.j]
   [1.25186920e+00+0.j]
   [1.25520528e+00+0.j]
   [1.25843261e+00+0.j]
   [1.26155796e+00+0.j]
   [1.26458744e+00+0.j]
   [1.26752657e+00+0.j]
   [1.27038039e+00+0.j]
   [1.27315350e+00+0.j]
   [1.27585010e+00+0.j]
   [1.27847405e+00+0.j]
   [1.28102889e+00+0.j]
   [1.28351789e+00+0.j]
   [1.28594406e+00+0.j]
   [1.28831021e+00+0.j]
   [1.29061893e+00+0.j]
   [1.29287263e+00+0.j]
   [1.29507354e+00+0.j]
   [1.29722377e+00+0.j]
   [1.29932526e+00+0.j]
   [1.30137985e+00+0.j]
   [1.30338924e+00+0.j]
   [1.30535505e+00+0.j]
   [1.30727878e+00+0.j]
   [1.30916185e+00+0.j]
   [1.31100560e+00+0.j]
   [1.31281129e+00+0.j]
   [1.31458012e+00+0.j]
   [1.31631320e+00+0.j]
   [1.31801160e+00+0.j]
   [1.31967633e+00+0.j]
   [1.32130834e+00+0.j]
   [1.32290854e+00+0.j]
   [1.32447779e+00+0.j]
   [1.32601691e+00+0.j]
   [1.32752666e+00+0.j]
   [1.32900780e+00+0.j]
   [1.33046103e+00+0.j]
   [1.33188702e+00+0.j]
   [1.33328641e+00+0.j]
   [1.33465981e+00+0.j]
   [1.33600780e+00+0.j]
   [1.33733096e+00+0.j]
   [1.33862980e+00+0.j]
   [1.33990484e+00+0.j]
   [1.34115656e+00+0.j]
   [1.34238545e+00+0.j]
   [1.34359194e+00+0.j]
   [1.34477647e+00+0.j]
   [1.34593945e+00+0.j]
   [1.34708128e+00+0.j]
   [1.34820234e+00+0.j]
   [1.34930298e+00+0.j]
   [1.35038358e+00+0.j]
   [1.35144446e+00+0.j]
   [1.35248594e+00+0.j]
   [1.35350835e+00+0.j]
   [1.35451198e+00+0.j]
   [1.35549712e+00+0.j]
   [1.35646406e+00+0.j]
   [1.35741305e+00+0.j]
   [1.35834437e+00+0.j]
   [1.35925825e+00+0.j]
   [1.36015495e+00+0.j]
   [1.36103469e+00+0.j]
   [1.36189770e+00+0.j]
   [1.36274419e+00+0.j]
   [1.36357437e+00+0.j]
   [1.36438845e+00+0.j]
   [1.36518661e+00+0.j]
   [1.36596906e+00+0.j]
   [1.36673596e+00+0.j]
   [1.36748750e+00+0.j]
   [1.36822384e+00+0.j]
   [1.36894514e+00+0.j]
   [1.36965157e+00+0.j]
   [1.37034328e+00+0.j]
   [1.37102042e+00+0.j]
   [1.37168312e+00+0.j]
   [1.37233153e+00+0.j]
   [1.37296577e+00+0.j]
   [1.37358599e+00+0.j]
   [1.37419230e+00+0.j]
   [1.37478482e+00+0.j]
   [1.37536368e+00+0.j]
   [1.37592897e+00+0.j]
   [1.37648081e+00+0.j]
   [1.37701932e+00+0.j]
   [1.37754457e+00+0.j]
   [1.37805668e+00+0.j]
   [1.37855575e+00+0.j]
   [1.37904185e+00+0.j]
   [1.37951508e+00+0.j]
   [1.37997552e+00+0.j]
   [1.38042326e+00+0.j]
   [1.38085838e+00+0.j]
   [1.38128094e+00+0.j]
   [1.38169103e+00+0.j]
   [1.38208872e+00+0.j]
   [1.38247406e+00+0.j]
   [1.38284714e+00+0.j]
   [1.38320801e+00+0.j]
   [1.38355672e+00+0.j]
   [1.38389335e+00+0.j]
   [1.38421795e+00+0.j]
   [1.38453056e+00+0.j]
   [1.38483124e+00+0.j]
   [1.38512004e+00+0.j]
   [1.38539701e+00+0.j]
   [1.38566219e+00+0.j]
   [1.38591561e+00+0.j]
   [1.38615734e+00+0.j]
   [1.38638739e+00+0.j]
   [1.38660581e+00+0.j]
   [1.38681263e+00+0.j]
   [1.38700788e+00+0.j]
   [1.38719160e+00+0.j]
   [1.38736382e+00+0.j]
   [1.38752455e+00+0.j]
   [1.38767383e+00+0.j]
   [1.38781167e+00+0.j]
   [1.38793810e+00+0.j]
   [1.38805314e+00+0.j]
   [1.38815681e+00+0.j]
   [1.38824912e+00+0.j]
   [1.38833008e+00+0.j]
   [1.38839971e+00+0.j]
   [1.38845802e+00+0.j]
   [1.38850502e+00+0.j]
   [1.38854071e+00+0.j]
   [1.38856510e+00+0.j]
   [1.38857819e+00+0.j]
   [1.38857999e+00+0.j]
   [1.38857049e+00+0.j]
   [1.38854970e+00+0.j]
   [1.38851761e+00+0.j]
   [1.38847422e+00+0.j]
   [1.38841951e+00+0.j]
   [1.38835349e+00+0.j]
   [1.38827614e+00+0.j]
   [1.38818745e+00+0.j]
   [1.38808740e+00+0.j]
   [1.38797598e+00+0.j]
   [1.38785318e+00+0.j]
   [1.38771897e+00+0.j]
   [1.38757334e+00+0.j]
   [1.38741626e+00+0.j]
   [1.38724771e+00+0.j]
   [1.38706766e+00+0.j]
   [1.38687608e+00+0.j]
   [1.38667294e+00+0.j]
   [1.38645822e+00+0.j]
   [1.38623188e+00+0.j]
   [1.38599388e+00+0.j]
   [1.38574418e+00+0.j]
   [1.38548275e+00+0.j]
   [1.38520955e+00+0.j]
   [1.38492452e+00+0.j]
   [1.38462763e+00+0.j]
   [1.38431882e+00+0.j]
   [1.38399805e+00+0.j]
   [1.38366526e+00+0.j]
   [1.38332040e+00+0.j]
   [1.38296341e+00+0.j]
   [1.38259423e+00+0.j]
   [1.38221280e+00+0.j]
   [1.38181905e+00+0.j]
   [1.38141292e+00+0.j]
   [1.38099434e+00+0.j]
   [1.38056323e+00+0.j]
   [1.38011952e+00+0.j]
   [1.37966313e+00+0.j]
   [1.37919398e+00+0.j]
   [1.37871199e+00+0.j]
   [1.37821706e+00+0.j]
   [1.37770912e+00+0.j]
   [1.37718806e+00+0.j]
   [1.37665379e+00+0.j]
   [1.37610620e+00+0.j]
   [1.37554521e+00+0.j]
   [1.37497068e+00+0.j]
   [1.37438253e+00+0.j]
   [1.37378062e+00+0.j]
   [1.37316485e+00+0.j]
   [1.37253508e+00+0.j]
   [1.37189120e+00+0.j]
   [1.37123307e+00+0.j]
   [1.37056054e+00+0.j]
   [1.36987349e+00+0.j]
   [1.36917177e+00+0.j]
   [1.36845522e+00+0.j]
   [1.36772368e+00+0.j]
   [1.36697701e+00+0.j]
   [1.36621502e+00+0.j]
   [1.36543754e+00+0.j]
   [1.36464441e+00+0.j]
   [1.36383542e+00+0.j]
   [1.36301039e+00+0.j]
   [1.36216911e+00+0.j]
   [1.36131139e+00+0.j]
   [1.36043700e+00+0.j]
   [1.35954573e+00+0.j]
   [1.35863735e+00+0.j]
   [1.35771161e+00+0.j]
   [1.35676828e+00+0.j]
   [1.35580709e+00+0.j]
   [1.35482777e+00+0.j]
   [1.35383006e+00+0.j]
   [1.35281367e+00+0.j]
   [1.35177829e+00+0.j]
   [1.35072362e+00+0.j]
   [1.34964935e+00+0.j]
   [1.34855512e+00+0.j]
   [1.34744061e+00+0.j]
   [1.34630544e+00+0.j]
   [1.34514924e+00+0.j]
   [1.34397162e+00+0.j]
   [1.34277216e+00+0.j]
   [1.34155046e+00+0.j]
   [1.34030606e+00+0.j]
   [1.33903849e+00+0.j]
   [1.33774729e+00+0.j]
   [1.33643194e+00+0.j]
   [1.33509192e+00+0.j]
   [1.33372667e+00+0.j]
   [1.33233562e+00+0.j]
   [1.33091818e+00+0.j]
   [1.32947370e+00+0.j]
   [1.32800152e+00+0.j]
   [1.32650096e+00+0.j]
   [1.32497128e+00+0.j]
   [1.32341172e+00+0.j]
   [1.32182147e+00+0.j]
   [1.32019968e+00+0.j]
   [1.31854548e+00+0.j]
   [1.31685791e+00+0.j]
   [1.31513599e+00+0.j]
   [1.31337867e+00+0.j]
   [1.31158484e+00+0.j]
   [1.30975335e+00+0.j]
   [1.30788294e+00+0.j]
   [1.30597231e+00+0.j]
   [1.30402007e+00+0.j]
   [1.30202473e+00+0.j]
   [1.29998472e+00+0.j]
   [1.29789836e+00+0.j]
   [1.29576386e+00+0.j]
   [1.29357931e+00+0.j]
   [1.29134267e+00+0.j]
   [1.28905173e+00+0.j]
   [1.28670414e+00+0.j]
   [1.28429738e+00+0.j]
   [1.28182871e+00+0.j]]
  </code></pre>
  <p>This is the current value of theta:</p>
  <pre><code>[-3.14159265e+00 -3.13159265e+00 -3.12159265e+00 -3.11159265e+00
   -3.10159265e+00 -3.09159265e+00 -3.08159265e+00 -3.07159265e+00
   -3.06159265e+00 -3.05159265e+00 -3.04159265e+00 -3.03159265e+00
   -3.02159265e+00 -3.01159265e+00 -3.00159265e+00 -2.99159265e+00
   -2.98159265e+00 -2.97159265e+00 -2.96159265e+00 -2.95159265e+00
   -2.94159265e+00 -2.93159265e+00 -2.92159265e+00 -2.91159265e+00
   -2.90159265e+00 -2.89159265e+00 -2.88159265e+00 -2.87159265e+00
   -2.86159265e+00 -2.85159265e+00 -2.84159265e+00 -2.83159265e+00
   -2.82159265e+00 -2.81159265e+00 -2.80159265e+00 -2.79159265e+00
   -2.78159265e+00 -2.77159265e+00 -2.76159265e+00 -2.75159265e+00
   -2.74159265e+00 -2.73159265e+00 -2.72159265e+00 -2.71159265e+00
   -2.70159265e+00 -2.69159265e+00 -2.68159265e+00 -2.67159265e+00
   -2.66159265e+00 -2.65159265e+00 -2.64159265e+00 -2.63159265e+00
   -2.62159265e+00 -2.61159265e+00 -2.60159265e+00 -2.59159265e+00
   -2.58159265e+00 -2.57159265e+00 -2.56159265e+00 -2.55159265e+00
   -2.54159265e+00 -2.53159265e+00 -2.52159265e+00 -2.51159265e+00
   -2.50159265e+00 -2.49159265e+00 -2.48159265e+00 -2.47159265e+00
   -2.46159265e+00 -2.45159265e+00 -2.44159265e+00 -2.43159265e+00
   -2.42159265e+00 -2.41159265e+00 -2.40159265e+00 -2.39159265e+00
   -2.38159265e+00 -2.37159265e+00 -2.36159265e+00 -2.35159265e+00
   -2.34159265e+00 -2.33159265e+00 -2.32159265e+00 -2.31159265e+00
   -2.30159265e+00 -2.29159265e+00 -2.28159265e+00 -2.27159265e+00
   -2.26159265e+00 -2.25159265e+00 -2.24159265e+00 -2.23159265e+00
   -2.22159265e+00 -2.21159265e+00 -2.20159265e+00 -2.19159265e+00
   -2.18159265e+00 -2.17159265e+00 -2.16159265e+00 -2.15159265e+00
   -2.14159265e+00 -2.13159265e+00 -2.12159265e+00 -2.11159265e+00
   -2.10159265e+00 -2.09159265e+00 -2.08159265e+00 -2.07159265e+00
   -2.06159265e+00 -2.05159265e+00 -2.04159265e+00 -2.03159265e+00
   -2.02159265e+00 -2.01159265e+00 -2.00159265e+00 -1.99159265e+00
   -1.98159265e+00 -1.97159265e+00 -1.96159265e+00 -1.95159265e+00
   -1.94159265e+00 -1.93159265e+00 -1.92159265e+00 -1.91159265e+00
   -1.90159265e+00 -1.89159265e+00 -1.88159265e+00 -1.87159265e+00
   -1.86159265e+00 -1.85159265e+00 -1.84159265e+00 -1.83159265e+00
   -1.82159265e+00 -1.81159265e+00 -1.80159265e+00 -1.79159265e+00
   -1.78159265e+00 -1.77159265e+00 -1.76159265e+00 -1.75159265e+00
   -1.74159265e+00 -1.73159265e+00 -1.72159265e+00 -1.71159265e+00
   -1.70159265e+00 -1.69159265e+00 -1.68159265e+00 -1.67159265e+00
   -1.66159265e+00 -1.65159265e+00 -1.64159265e+00 -1.63159265e+00
   -1.62159265e+00 -1.61159265e+00 -1.60159265e+00 -1.59159265e+00
   -1.58159265e+00 -1.57159265e+00 -1.56159265e+00 -1.55159265e+00
   -1.54159265e+00 -1.53159265e+00 -1.52159265e+00 -1.51159265e+00
   -1.50159265e+00 -1.49159265e+00 -1.48159265e+00 -1.47159265e+00
   -1.46159265e+00 -1.45159265e+00 -1.44159265e+00 -1.43159265e+00
   -1.42159265e+00 -1.41159265e+00 -1.40159265e+00 -1.39159265e+00
   -1.38159265e+00 -1.37159265e+00 -1.36159265e+00 -1.35159265e+00
   -1.34159265e+00 -1.33159265e+00 -1.32159265e+00 -1.31159265e+00
   -1.30159265e+00 -1.29159265e+00 -1.28159265e+00 -1.27159265e+00
   -1.26159265e+00 -1.25159265e+00 -1.24159265e+00 -1.23159265e+00
   -1.22159265e+00 -1.21159265e+00 -1.20159265e+00 -1.19159265e+00
   -1.18159265e+00 -1.17159265e+00 -1.16159265e+00 -1.15159265e+00
   -1.14159265e+00 -1.13159265e+00 -1.12159265e+00 -1.11159265e+00
   -1.10159265e+00 -1.09159265e+00 -1.08159265e+00 -1.07159265e+00
   -1.06159265e+00 -1.05159265e+00 -1.04159265e+00 -1.03159265e+00
   -1.02159265e+00 -1.01159265e+00 -1.00159265e+00 -9.91592654e-01
   -9.81592654e-01 -9.71592654e-01 -9.61592654e-01 -9.51592654e-01
   -9.41592654e-01 -9.31592654e-01 -9.21592654e-01 -9.11592654e-01
   -9.01592654e-01 -8.91592654e-01 -8.81592654e-01 -8.71592654e-01
   -8.61592654e-01 -8.51592654e-01 -8.41592654e-01 -8.31592654e-01
   -8.21592654e-01 -8.11592654e-01 -8.01592654e-01 -7.91592654e-01
   -7.81592654e-01 -7.71592654e-01 -7.61592654e-01 -7.51592654e-01
   -7.41592654e-01 -7.31592654e-01 -7.21592654e-01 -7.11592654e-01
   -7.01592654e-01 -6.91592654e-01 -6.81592654e-01 -6.71592654e-01
   -6.61592654e-01 -6.51592654e-01 -6.41592654e-01 -6.31592654e-01
   -6.21592654e-01 -6.11592654e-01 -6.01592654e-01 -5.91592654e-01
   -5.81592654e-01 -5.71592654e-01 -5.61592654e-01 -5.51592654e-01
   -5.41592654e-01 -5.31592654e-01 -5.21592654e-01 -5.11592654e-01
   -5.01592654e-01 -4.91592654e-01 -4.81592654e-01 -4.71592654e-01
   -4.61592654e-01 -4.51592654e-01 -4.41592654e-01 -4.31592654e-01
   -4.21592654e-01 -4.11592654e-01 -4.01592654e-01 -3.91592654e-01
   -3.81592654e-01 -3.71592654e-01 -3.61592654e-01 -3.51592654e-01
   -3.41592654e-01 -3.31592654e-01 -3.21592654e-01 -3.11592654e-01
   -3.01592654e-01 -2.91592654e-01 -2.81592654e-01 -2.71592654e-01
   -2.61592654e-01 -2.51592654e-01 -2.41592654e-01 -2.31592654e-01
   -2.21592654e-01 -2.11592654e-01 -2.01592654e-01 -1.91592654e-01
   -1.81592654e-01 -1.71592654e-01 -1.61592654e-01 -1.51592654e-01
   -1.41592654e-01 -1.31592654e-01 -1.21592654e-01 -1.11592654e-01
   -1.01592654e-01 -9.15926536e-02 -8.15926536e-02 -7.15926536e-02
   -6.15926536e-02 -5.15926536e-02 -4.15926536e-02 -3.15926536e-02
   -2.15926536e-02 -1.15926536e-02 -1.59265359e-03  8.40734641e-03
    1.84073464e-02  2.84073464e-02  3.84073464e-02  4.84073464e-02
    5.84073464e-02  6.84073464e-02  7.84073464e-02  8.84073464e-02
    9.84073464e-02  1.08407346e-01  1.18407346e-01  1.28407346e-01
    1.38407346e-01  1.48407346e-01  1.58407346e-01  1.68407346e-01
    1.78407346e-01  1.88407346e-01  1.98407346e-01  2.08407346e-01
    2.18407346e-01  2.28407346e-01  2.38407346e-01  2.48407346e-01
    2.58407346e-01  2.68407346e-01  2.78407346e-01  2.88407346e-01
    2.98407346e-01  3.08407346e-01  3.18407346e-01  3.28407346e-01
    3.38407346e-01  3.48407346e-01  3.58407346e-01  3.68407346e-01
    3.78407346e-01  3.88407346e-01  3.98407346e-01  4.08407346e-01
    4.18407346e-01  4.28407346e-01  4.38407346e-01  4.48407346e-01
    4.58407346e-01  4.68407346e-01  4.78407346e-01  4.88407346e-01
    4.98407346e-01  5.08407346e-01  5.18407346e-01  5.28407346e-01
    5.38407346e-01  5.48407346e-01  5.58407346e-01  5.68407346e-01
    5.78407346e-01  5.88407346e-01  5.98407346e-01  6.08407346e-01
    6.18407346e-01  6.28407346e-01  6.38407346e-01  6.48407346e-01
    6.58407346e-01  6.68407346e-01  6.78407346e-01  6.88407346e-01
    6.98407346e-01  7.08407346e-01  7.18407346e-01  7.28407346e-01
    7.38407346e-01  7.48407346e-01  7.58407346e-01  7.68407346e-01
    7.78407346e-01  7.88407346e-01  7.98407346e-01  8.08407346e-01
    8.18407346e-01  8.28407346e-01  8.38407346e-01  8.48407346e-01
    8.58407346e-01  8.68407346e-01  8.78407346e-01  8.88407346e-01
    8.98407346e-01  9.08407346e-01  9.18407346e-01  9.28407346e-01
    9.38407346e-01  9.48407346e-01  9.58407346e-01  9.68407346e-01
    9.78407346e-01  9.88407346e-01  9.98407346e-01  1.00840735e+00
    1.01840735e+00  1.02840735e+00  1.03840735e+00  1.04840735e+00
    1.05840735e+00  1.06840735e+00  1.07840735e+00  1.08840735e+00
    1.09840735e+00  1.10840735e+00  1.11840735e+00  1.12840735e+00
    1.13840735e+00  1.14840735e+00  1.15840735e+00  1.16840735e+00
    1.17840735e+00  1.18840735e+00  1.19840735e+00]
  
  </code></pre>
  <p>The Matlab script converges in a splut of a second, while my algorithm seems to run for a very long time. The output error is very large since the value of r is too small (see output above).</p>