Ordinary Least Squares Regression with binary dependent variable I know that OLS regression is linear and output expected is continuous and values will fall higher than 1 or less than 0 so is no meaning of values what are not between 0 and 1 (here pointing to values 20, 30 etc not strictly around 1), can’t be interpreted.
My questions is: Mathematically speaking: Why is not appropriate to use OLS regression when we have binary dependent variable? 
 A: This approach can (and likely will) produce values that are impossible, and there are other methods like logistic regression that avoid this issue. That’s one reason why OLS might not be preferred for such a situation.
However, the math works out fine. Nothing will stop you from fitting an OLS model with your parameters estimated the usual way with $\hat\beta=(X^TX)^{-1}X^Ty$. If such a model works best for what you’re doing, as proponents of linear probability models believe, then go for it! Nothing about the derivation of $\hat\beta=(X^TX)^{-1}X^Ty$ makes a distribution assumption about the error term.
For instance, in the following simulation, the linear probability model has tight confidence intervals around the true parameters (intercept is $0$, slope is $1$), and that even happens with a small sample size of just ten.
set.seed(2023)
R <- 1000
N <- 10
intercepts <- slopes <- rep(NA, R)
for (i in 1:R){
  x <- runif(N, 0, 1)
  y <- rbinom(N, 1, x)
  L <- lm(y ~ x)
  intercepts[i] <- summary(L)$coef[1, 1]
  slopes[i] <- summary(L)$coef[2, 1]
}
t.test(intercepts)$conf.int[c(1, 2)] # Confidence interval (-0.01600467,  0.01664349)
t.test(slopes, mu = 1)$conf.int[c(1, 2)] # (0.9723355, 1.0277178)

If the true relationship is $\mathbb E\left[Y\vert X\right] = \beta_0 + \beta_1X$, as is the case here, then the linear model captures this relationship.
As with Frank Harrell in the comments, I have my doubts about how well you will get a linear probability model to fit most data, but if it works...
