Q:

How to run a regression for a given value of a parameter and represent this in a plot and stagazer?

Given the linear model:
    data(wage1)
    reg_wage1 <- lm(lwage ~ female + educ + (female * educ), data = wage1)

stargazer(reg_wage1, type = "text")
===============================================
                        Dependent variable:    
                    ---------------------------
                               lwage           
-----------------------------------------------
female                        -0.360*          
                              (0.185)          
                                               
educ                         0.077***          
                              (0.009)          
                                               
female:educ                   -0.0001          
                              (0.015)          
                                               
Constant                     0.826***          
                              (0.118)          
                                               
-----------------------------------------------
Observations                    526            
R2                             0.300           
Adjusted R2                    0.296           
Residual Std. Error      0.446 (df = 522)      
F Statistic           74.649*** (df = 3; 522)  
===============================================
Note:               *p<0.1; **p<0.05; ***p<0.01

Create a line graph that compares each gender's wage growth (female is a data type that can only be 1 or 0. 1 means yes this is a female and 0 means no.), illustrates two graphs corresponding to different values of $delta_0$ and $delta_1$.
My attempt:
    plot(NULL, xlim=c(0,16), ylim=c(-1,1), ylab="Log Wage", xlab="Education", main = "Men vs Women's Wages Slope")
    abline(a = -0.360, b = .185, col = "blue")
    abline(a = 0.077, b = 0.009, col = "red")
    legend("bottomright", c("Female", "Education"),
             col = c("blue", "red"), fill= c("blue", "red"), title="Legend", cex=.75)

Where a(constant) is the intercept and b(standard error) is the slope (these values are taken from the stargazer table).
My main problem is I don't know how to run the model where I can evaluate female and male wage growth independently, that is where female = 1 or female = 0. I tried using predict() but the result it returns does not come in a format where I am able to get the coefficients in I need to plot into plot()
The results of this exercise should look something like this:

A:

First of all, please make sure to note where your data comes from (in this case the np package).
Your intuition of using the predict function was one way of doing it but I would recommend using the ggplot2 package to plot your regression results. 
The easiest implementation of this is using the geom_smooth argument, which does the modelling automatically for you and because we specify a color argument, your interaction is taken into account: 
library(np)
data(wage1)

ggplot(wage1, aes(x=educ, y=lwage, color=female))+
  geom_point()+
  geom_smooth(method = "lm")

You could also use the predict function, as mentioned, but you have to run the model manually before (and adding the confidence intervals is a bit of a pain):
reg_wage1 <- lm(lwage ~ female + educ + (female * educ), data = wage1)
wage1$pred <- predict(reg_wage1)

ggplot(wage1, aes(x=educ, y=lwage, color=female))+
  geom_point()+
  geom_line(aes(x = educ, y = pred, color=female))

Both lines have different slopes as well, although they are really difficult to see with your model. To reassure you that this is not just an intercept correction but also shows you different slopes, consider another example where the interaction effects are much stronger: 
data(mtcars)
mtcars$vs <- as.factor(mtcars$vs)
ggplot(mtcars,aes(x=hp,y=mpg,color=vs)) + 
  geom_point() + 
  geom_smooth(method = "lm")

Hope this helps!

