Conceptually, $\epsilon_i$ is the difference between an individual $i$'s actual wage and your model's prediction for $i$'s wage.
For example, let's assume that in the data, individual $i$ has $e_i = 12$ years of education and a wage of $w_i = 32{,}250$. Then the model prediction for $i$'s wage (given the coefficients you listed) is $13{,}000 + 1{,}500*12 = 31{,}000$. But we observed that $i$'s actual wage is $32{,}250$. Then $\epsilon_i = 32{,}250 - 31{,}000 = 1{,}250$. This is just a concrete example of the formula in @roundsquare's comment.
The $\epsilon_i$ are illustrated with fake data in the plot below. (The R code to create the plot is just below the plot.) We fit the model $y_i = \beta_0 + \beta_1 x_i + \epsilon_i$. In the plot, the points are the data, the line is the regression fit (where the function is $y = \beta_0 + \beta_1 x$) and the red vertical lines are the $\epsilon_i$ (the difference between the regression line and the actual data), which are known as the residuals. The regression line is the particular line that minimizes $\Sigma{\epsilon_i^2}$ (the sum of the squared residuals).
library(ggplot2)
# Fake data
set.seed(491)
x = runif(50, 0, 100)
y = 2*x + 5 + rnorm(50,0,20)
# Combine x, y into a data frame
dat = data.frame(x,y)
# Fit regression model
m1 = lm(y ~ x, data=dat)
# Get predictions for regression model
dat$ypred = predict(m1)
ggplot(dat, aes(x,y)) +
geom_smooth(method="lm", se=FALSE, size=0.7) +
geom_segment(aes(yend=ypred, xend=x), colour="red") +
geom_point(pch=21, fill=hcl(240,100,55)) +
theme_bw()
In R, you can extract all of the $\epsilon_i$ from the regression model m1 above by running resid(m1) or dat$y - predict(m1).
