mcoxph<-function(formula){
	 hh<-0.0001
	#### we use the same function arguments as coxph, so we can call coxph
   ### steal first few lines from glm
	m <- match.call(expand = F)
	m[[1]] <- as.name("coxph")
	coxobj <- eval(m, sys.parent(1))
	#coxobj<-coxph(Surv(time,status)~x)
beta<-coxobj$coefficients
scorei<-residuals.coxph(coxobj,type="score")
N<-coxobj$n
nbeta<- length(as.vector(beta))
scorep<-array(0,c(N,nbeta,nbeta))
scorem<-array(0,c(N,nbeta,nbeta))
zeros<-rep(0,nbeta)
for (i in 1:nbeta){
	delta<- zeros
	delta[i]<- hh
		callval<- call("coxph",formula=m$formula,init=beta+delta,iter.max=0)
	coxobj <- eval(callval, sys.parent(2))
#   coxobj<-coxph(Surv(time,status)~x,init=beta+delta,iter.max=0)
   temp<-residuals.coxph(coxobj,type="score")
	scorep[,i,]<- scorep[,i,] + temp
	scorep[,,i]<-scorep[,,i] + temp
	callval<- call("coxph",formula=m$formula,init=beta-delta,iter.max=0)
	coxobj <- eval(callval, sys.parent(2))
	#coxobj<-coxph(Surv(time,status)~x,init=beta-delta,iter.max=0)

   temp<-residuals.coxph(coxobj,type="score")
	scorem[,i,]<- scorem[,i,] + temp
	scorem[,,i]<-scorem[,,i] + temp
	}

omega<- -(scorep - scorem)/(4*hh)
list(beta=beta,u=scorei,omega=omega)
}