function P = SmoothSplineTol(G,tol)

beta = 1000;
epsilon = 10^(-4);
[P,C,A] = SmoothSpline2(G,beta);
M = A'*A;
Fs = A\G'-P;
d = Fs'*M*Fs;

while (abs(tol-d) > epsilon)  
    dd = -2*Fs'*M*(C+beta*M)^(-1)*M*Fs;
    q = -d/dd-beta;
    p = -d^2/dd;
    beta = (p-q*tol)/tol;
    [P,C,A] = SmoothSpline2(G,beta);
    Fs = A\G'-P;
    d = Fs'*M*Fs;
    tol-d
end

beta