ml command error "could not find feasible values" with a seemingly non-problematic mata-based evaluator

[Minchul Park]

I am writing my own tobit5 program. But, the ml command returns the error "could not find feasible values".
Specifically, my do file (not an ado file) is like below:

program define my_tobit5, eclass
    version 18
    
    syntax varlist(min = 1) [if] [in], D(varname)
    
    gettoken y xvar : varlist
    gettoken D : d
    
    ml model lf my_llk() (`D' = `xvar') (`y' = `xvar') (`y' = `xvar') () () () ()
    ml check
    ml maximize
end

version 18
mata:
mata set matastrict on
void my_llk(transmorphic scalar ML,
             real rowvector b,
             real colvector lnf)
{
    real vector D, depvar, xgam, xb0, xb1, u0, u1, llk0, llk1, ///
    sig0, sig1, rho0, rho1, rho0sq, rho1sq
    
    D      = moptimize_util_depvar(ML, 1)
    xgam   = moptimize_util_xb(ML, b, 1)
    
    depvar = moptimize_util_depvar(ML, 2)
    xb0    = moptimize_util_xb(ML, b, 2)
    xb1    = moptimize_util_xb(ML, b, 3)
    sig0   = moptimize_util_xb(ML, b, 4)
    sig1   = moptimize_util_xb(ML, b, 5)
    rho0   = moptimize_util_xb(ML, b, 6)
    rho1   = moptimize_util_xb(ML, b, 7)
        
    u0     = depvar :- xb0
    u1     = depvar :- xb1
    rho0sq = sqrt(1 :- rho0:^2)
    rho1sq = sqrt(1 :- rho1:^2)
    
    llk0   = lnnormal((-xgam :- (rho0 :/ sig0) :* u0) :/ rho0sq) :- ln(sig0) :- 0.5 :* (u0 :/ sig0):^2
    llk1   = lnnormal((xgam :+ (rho1 :/ sig1) :* u1) :/ rho1sq) :- ln(sig1) :- 0.5 :* (u1 :/ sig1):^2
    
    lnf    = ((1 :- D) :* llk0) :+ (D :* llk1)
}
end

That is, the first program my_tobit5 includes the ml command that evaluate a mata-based evaluator my_llk.
When I run the do file and use the command below, I get the error "could not find feasible values".

sysuse auto
my_tobit5 price mpg, d(foreign)

Since the return shows that the first evaluated value of the likelihood function is minus infinite and there is no other feasible value,
I checked whether the mata function my_llk returns minus infinite always by using a pseudo data sample.
But, the return lnf is reasonably just a vector of real numbers.

So, what is the problem?

p.s. the following example shows that our mata function returns individual likelihood function values.
(the data generating process is not consistent with the tobit model, it is just a example DGP)

mata:
D      = J(250, 1, (0 \ 1))
X      = (J(500, 1, 1), rnormal(500, 3, 0, 1))
gam    = (2 \ 1 \ -1 \ 2)
b0     = (-1 \ 1 \ 1 \ -1)
b1     = (1 \ 0 \ 1 \ 1)
xgam   = X * gam
xb0    = X * b0
xb1    = X * b1
depvar = xb1 + rnormal(500, 1, 0, 1)

sig0   = J(500, 1, 0.5)
sig1   = J(500, 1, 0.2)
rho0   = J(500, 1, -0.2)
rho1   = J(500, 1, 0.2)

u0     = depvar :- xb0
u1     = depvar :- xb1
rho0sq = sqrt(1 :- rho0:^2)
rho1sq = sqrt(1 :- rho1:^2)

llk0   = lnnormal((-xgam :- (rho0 :/ sig0) :* u0) :/ rho0sq) :- ln(sig0) :- 0.5 :* (u0 :/ sig0):^2
llk1   = lnnormal((xgam :+ (rho1 :/ sig1) :* u1) :/ rho1sq) :- ln(sig1) :- 0.5 :* (u1 :/ sig1):^2

lnf    = ((1 :- D) :* llk0) :+ (D :* llk1)

lnf
end