In general you would need a solver. One exception would be a fully saturated model, that is, a model with only categorical explanatory/right-hand-side/x-variables and one parameter for each combination of groups. In that case the coefficients are just a function of the group means. In case of a bivariate model, it would mean that your only explanatory/right-hand-side/x-variable is a binary variable. Below is an example of how to recreate the results obtained by a solver of a fully saturated bivariate model by just transforming means in Stata.
// load some example data
sysuse nlsw88, clear
// use a "solver"
logit union collgrad
// collect the means
tempname noncoll coll
sum union if e(sample) & collgrad == 0, meanonly
scalar `noncoll' = r(mean)
sum union if e(sample) & collgrad == 1, meanonly
scalar `coll' = r(mean)
// recreate the coefficients using just means
di as txt "the constant" ///
as result logit(`noncoll')
di as txt "coefficient of collgrad: " ///
as result logit(`coll') - logit(`noncoll')