%%%%%%%%%%% maxcurve %%%%%%%%%%
%
%find maximum curvature point in L-curve for Tikhonov Regularization

function out = maxcurve(tau,Lambda,u)
temp1 = diag(Lambda)+tau;
temp2 = 1./(temp1).^2;
temp3 = 1./(temp1).^3;
temp4 = 1./(temp1).^4;
u2 = u.*u;
out = -6*(sum(temp4.*u2))*(1+(sum(temp2.*u2))^2) + ...
    12*(sum(temp3.*u2))^2*(sum(temp2.*u2));
end