your code as provided errors, there are undeclared objects such as aa, bb etc.
assuming you have working code, the result of the code is a data.frame
If you want to save the data.frame to disk you can use saveRDS
whose parameters are first object, then file (as in filename).
i.e.

I have a long code. This code not the real code, but I want like the following
s1 s2 The True P-value F-SPA S-SPA T-SPA
1 3 --------- ---- ----- ------
1 4 --------- ---- ----- ------
2 3 --------- ---- ----- ------
2 4 --------- ---- ----- ------

How I can save like this in dataframe and excel file.

I don't understand what you are asking.
Does the code you have work, and you want to save the result (which I explained how to do)
or does your code not even work?

# clear data frame
rm(list = ls())
#-------------------------------------------------------
options(scipen=999) # turn-off scientific notation like 1e+48
#------------------------------------------------------------
lambda1=2
lambda2=1
lambda3=1
countt=0
ll=0
for (s1 in 1:10)
{
for (s2 in 1:10)
{
# The exact probability
for (j in 1:10000000)
{
f = rpois(1,lambda1)
if(f == 0)
{
countt = countt + 1
next
}
xx = rpois(f,lambda2)
yy = rpois(f,lambda3)
if(sum(xx)== s1 && sum(yy)<= s2)
{
ll = ll + 1
}
if(sum(xx)< s1 && sum(yy)== s2)
{
ll = ll + 1
}
if(sum(xx)< s1 && sum(yy)< s2)
{
countt=countt+1
}
}
pp1=(countt)/10000000
pp2=(countt+ll)/10000000
pp =(countt+0.5*ll)/10000000
# The saddlepoint Approximation
a=(s1+s2)*exp(lambda2+lambda3)/lambda1
t=log(s1*lambertW0(a)/(lambda2*(s1+s2)))
u=log(s2*lambertW0(a)/(lambda3*(s1+s2)))
b=s1*exp(lambda2)/lambda1
th=log(lambertW0(b)/(lambda2))
k1_th = lambda1*exp(lambda2*(exp(th)-1))*exp(lambda3*(exp(0)-1))-lambda1
T1 = sign(th)*sqrt(-2*(-th*s1+k1_th))
ktu = lambda1*exp(lambda2*(exp(t)-1))*exp(lambda3*(exp(u)-1))-lambda1
v0 = sign(u)*sqrt(-2*(ktu-k1_th-(t-th)*s1-u*s2))
k2_sh = lambda1*exp(lambda2*(exp(0)-1))*exp(lambda3*(exp(u)-1))-lambda1
w0 = sign(t)*sqrt(-2*(-t*s1-k2_sh+ktu))
b = (w0-T1)/v0
r1 = -b/sqrt(b^2+1)
v1 = (-b*T1+v0)/sqrt(b^2+1)
kuu = lambda1*exp(lambda2*(exp(t)-1))*lambda3*exp(u)*exp(lambda3*(exp(u)-1))+lambda1*exp(lambda2*(exp(t)-1))*lambda3^2*(exp(u))^2*exp(lambda3*(exp(u)-1))
ktt = lambda1*lambda2*exp(t)*exp(lambda2*(exp(t)-1))*exp(lambda3*(exp(u)-1))+lambda1*lambda2^2*(exp(t))^2*exp(lambda2*(exp(t)-1))*exp(lambda3*(exp(u)-1))
ksstt = lambda1*lambda2*exp(t)*exp(lambda2*(exp(t)-1))*lambda3*exp(u)*exp(lambda3*(exp(u)-1))
G = sqrt(kuu-ksstt^2/ktt)
I11=pbivnorm(T1,v1,r1)
I12=pnorm(w0)*dnorm(v0)*(1/v0-(u*G)^(-1))
I21=pnorm(v0)*dnorm(T1)*(1/w0-1/(t*ktt^(1/2)))
I22= exp(-t*s1-u*s2+ktu)*(1/w0-1/(t*ktt^(1/2)))*(1/v0-(u*G)^(-1))/(2*pi)
F1= I11+I12+I21+I22
I21d = pnorm(v0)*dnorm(T1)*(1/w0-1/((1-exp(-t))*ktt^(1/2)))
I22d = exp(-t*s1-u*s2+ktu)*(1/w0-1/((1-exp(-t))*ktt^(1/2)))*(1/v0-(u*G)^(-1))/(2*pi)
F2 = I11+I12+I21d+I22d
I12d = pnorm(w0)*dnorm(v0)*(1/v0-((1-exp(-u))*G)^(-1))
I22dd = exp(-t*s1-u*s2+ktu)*(1/w0-1/((1-exp(-t))*ktt^(1/2)))*(1/v0-(u*G)^(-1))/(2*pi)
F3 = I11+I12d+I21d+I22dd
print(paste0("s1 =",s1))
print(paste0("s2 =",s2))
print(paste0("The true p-value =",pp2))
print(paste0("The First approximation of Saddlepoint =",F1))
print(paste0("The Second approximation of Saddlepoint =",F2))
print(paste0("The Third approximation of Saddlepoint =",F3))
countt=0
ll=0
}
}

My code works. But I don't know how I save the result in dataframe or excel file. Can you help my?
I want to save the following
print(paste0("s1 =",s1))
print(paste0("s2 =",s2))
print(paste0("The true p-value =",pp2))
print(paste0("The F-SPA =",F1))
print(paste0("The S-SPA =",F2))
print(paste0("The T-SPA=",F3))

# I think this is what you have
for (i in 1:2){
print(paste0("i = ", i))
print(paste0("y = ", i*i))
}
# my suggestion
(gather_df <- data.frame(i=numeric(0),
y=numeric(0)))
for (i in 1:2){
gather_df <- rbind(gather_df,
data.frame(i=i,
y=i*i))
}
gather_df