secular<-function(b,a,k=2,xmin=-max(a),xmax=max(a),yup=100,ylw=0,off=1,eps=1e-6,itmax=100,verbose=FALSE) {
lb<-sqrt(2/k)*a; itel=1; xold<-off+max(0,max(-lb))
secPlot(b,lb,xmin,xmax,ylw,yup)
title(main=paste(deparse(b),deparse(a),sep=" "))
repeat {
	ff<-sum((b^2)/((lb+xold)^2))-(xold^2)
	gg<--2*(sum((b^2)/((lb+xold)^3))+xold)
	xnew<-xold-(ff/gg)
	if (verbose) cat(
		"Function: ",formatC(ff,digits=6,width=12,format="f"),
		"Gradient: ",formatC(gg,digits=6,width=12,format="f"),
		"Oldval:   ",formatC(xold,digits=6,width=12,format="f"),
		"Newval:   ",formatC(xnew,digits=6,width=12,format="f"),
		"\n")
	if ((abs(ff) < eps) || (itel == itmax)) break
	xold<-xnew; itel<-itel+1
	}
tht<-xnew; lbd<--tht^2; tau<-sqrt(2/k)*tht
y<--b/(tht*(lb+tht)); x<-sqrt(2/k)*tht*y
alpha<-sum(a*y*y); beta<-sum(b*y); bl<-beta+(tau*alpha)
cat("theta:  ",formatC(tht,digits=6,width=12,format="f"),"\n")
cat("lambda: ",formatC(lbd,digits=6,width=12,format="f"),"\n")
cat("tau:    ",formatC(tau,digits=6,width=12,format="f"),"\n")
cat("ssq:    ",formatC(sum(y^2),digits=6,width=12,format="f"),"\n")
cat("beta:   ",formatC(beta,digits=6,width=12,format="f"),"\n")
cat("ta+b:   ",formatC(bl,digits=6,width=12,format="f"),"\n")
cat("quadr:  ",formatC(bl+(k*tau*tau*.5),digits=6,width=12,format="f"),"\n")
cat("func:   ",formatC((tau*beta)+(alpha*tau*tau*.5)+(k*tau*tau*tau/6),digits=6,width=12,format="f"),"\n")
}

secPlot<-function(b,a,xmin,xmax,ylw,yup) {
x<-seq(xmin,xmax,length=1000)
y<-sapply(x,function(h) sum((b/(a+h))^2))
plot(x,y,ylim=c(ylw,yup))
points(x,x^2,col="RED")
}