ガウスニュートン法による関数フィッティング

function [serror,yfit]=gaussnewton(xdata,ydata,error)
syms a b x y

f = a*x/(x+b);
R = y - f;
S = R^2;

Jacob=subs(jacobian(R,[a,b]),{x,y},{xdata,ydata}); 
v = [0.9;0.2];
iteration = 100;
SEF = subs(S,{x,y},{xdata,ydata});

r = subs(R,{x,y},{xdata,ydata});

for i = 1:iteration
   J = double(subs(Jacob,[a,b],v'));
   df = double(subs(r,[a,b],v'));
   delta = - ( J.' * J )\J.' *  df;
   v = v + delta;

   SOE = double(subs(SEF,[a,b],v'));
   temp = SOE.' * SOE;
   if i==1
       serror=temp;
   else
       serror=horzcat(serror,temp);
   end
   clear temp
   if error>=serror(i)
       break;
   end
end

yfit = double( subs(subs(f,x,xdata),[a,b],v') );

end