Calculate best compounding period (cryptocurrency)

I have an interesting problem, I am not sure if the best place to put this in here or in stack mathematics but since I tried to solve it programmatically (in R, but any programming language would work) maybe StackOverflow is the best place. The problem involves compounding cryptocurrencies. The problem is the following: You are staking a "K" amount of token in a liquidity pool to get interests. The liquidity pool gives an "Annual Percentage Rate" (APR) to your token, i.e. an annual interest, the APR is not compounded. You can compound your interests any time you want but, every time you do it you have to pay a small fee.

I initially tried to solve it with some for loops, estimating the final return if a hypothetical user was compounding every "D" days. The current example uses CAKE and BNB tokens. A few assumptions to simplify the solution. Imagine that the APR, the bnb_value, and the cake_value are FIXED values (they are not in reality).

APR=1.349  ## the interest value showed in the pool, divided by 100
APR_day=APR/365.0 ## the daily APR
bnb_value=210 ## current value of the BNB token, in euro or dollar or any FIAT
fee_bnb=0.002*bnb_value ## current fee in euro or dollar or any FIAT
initial_cakes=10
cake_value=10

adf=data.frame(col1=NULL, col2=NULL, col3=NULL, col4=NULL)
## we generate a sequence containing fractional days up to the third day
## and then full days starting from the fourth day
comp_intervals=c(seq(0.1, 3, 0.1), 4:30)
for(j in comp_intervals){
  acquired_int=0  ## acquired interest in "j" days
  current_val=initial_cakes  ## current value of the capital
  all_fees=0  ## fees paid for the transactions
  aseq=seq(j, 365, j)  ## list of the compounding events in days
  
  # apr_seq is APR for the "j" period. If "j" is 1 it's the daily APR,
  # but for longer compounding intervals the APR for the "j" period
  # become larger
  apr_seq=APR_day*j
  for(i in aseq){
    acquired_int=current_val*apr_seq
    current_val=current_val+acquired_int
    all_fees=all_fees+fee_bnb
    acquired_int=0
  }
  
  ## we add the interest for the remaining days of the year, if any, to current_val
  acquired_int=(365-max(aseq))*APR_day*current_val
  current_val=current_val+acquired_int
  
  final_gain=round(current_val*cake_value - all_fees, 2)
  # msg=paste0("Final gain in Euro minus fees: ", final_gain)
  # print(msg)
  apy_i_day= round(final_gain/(initial_cakes*cake_value), 5)
  # msg=paste0("apy compounding every ", j, " days is: ", apy_i_day)
  # print(msg)
  # cat("\n")
  adf=rbind(adf, data.frame(fiat_value=final_gain,
      APY_val=100*apy_i_day, compounding_days=j, cakes=current_val))
  
}
# finally we show, who, among the various compounding, had the hiest yield
adf[adf$APY_val==max(adf$APY_val), ]

The problem is that the code you just saw doesn't really tell you what is the best period to compound your interests. It tells more what are the yields if a user compound every "D" days. It's close to the real solution, but it's not it! You can intuitively understand that it's wrong extending the time. You start with a small capital, so to get a good return you compound "rarely" because of the fees, but more time passes and more your capital grows. The more your capital grows more often you should compound.

I tried a different approach then.

The general formula that given the "Annual percentage yield" is the following:

APY = (1 + APR/N)^N -1

if you consider also the fees and the initial capital you have:

Final_capital=Initial_capital*APY - single_fee*N

Where APR is the Percentage rate, N is the number of compounding events (in this formula they are temporally equally distributed).

By differentiating Final_capital by "dN" and finding the zeros of the equation you get the best number of compounding events. If you divide 365 by the best number of compounding events you should get after how many days you should compound your tokens. The results I obtain from the differential formula are different from the first solution, I am not sure why. I also think, but I am not sure, that this latter solution has the same limitation of the previous one.

library(utils)
### interest APY formula minus fees
final_v=function(x, APR, fee_bnb, initial_value){
  return(  initial_value*( (1+APR/x)^x - 1 ) -fee_bnb*x   ) 
}

## differential respect to X of the interest APY formula minus the fees
a_diff_comp=function(x, APR, fee_bnb, initial_value){
  return( (initial_value*( (APR/x + 1)^x )*(  log(APR/x + 1) - APR/(  ((1/x)+1)*x  ) ) ) - fee_bnb )
}

x=3:40
y=sapply(x, a_diff_comp, APR=APR, fee_bnb=fee_bnb, initial_value=(initial_cakes*cake_value))
plot(x,y)

y2=sapply(x, final_v, APR=APR, fee_bnb=fee_bnb, initial_value=(initial_cakes*cake_value))
plot(x,y2)

xmin <- uniroot(a_diff_comp, c(1, 100), tol = 0.000001, APR=APR, fee_bnb=fee_bnb, initial_value=(initial_cakes*cake_value))
xmin$root

So, how to calculate properly the best compounding interval?