Intersect points and polygons

I'm afraid you're making this more complicated for yourself. I have created reproducible code that illustrates what you should do step by step. See below:

EDIT: After seeing what you were really looking for, which is essentially the answer noted in the other stack exchange post you linked here, I have added a reproducible example that should be able to guide you. See below. Note, I am leaving my previous answer below this one just in case it is helpful.

library(sf)
library(tmap)
library(dplyr)

data(World) # free dataset that will load
countries = World %>%
  filter(name == "Afghanistan" | 
           name == "Turkmenistan" |
           name == "Iran") #Picking three countries in the world

pnts = data.frame(name=c("Point1","Point2","Point3","Point4"),latitude = c(34,35,36,37),longitude=c(63,65,62.5,66.5)) #Picking random points

sp_points = st_as_sf(pnts,coords = c('longitude',"latitude"))#make points spatial
st_crs(sp_points)= 4326 # Give the points a coordinate reference system (CRS)
sp_points=st_transform(sp_points,crs = st_crs(countries)) # Match the point and polygon CRS

tm_shape(countries)+
  tm_borders() +
  tm_text('name')+
  tm_shape(sp_points)+
  tm_dots(col='red',size=1.2)+
  tm_text('name',ymod = 1) #Creates the map below

sp_points$country_name <- apply(st_intersects(countries, sp_points, sparse = FALSE), 2, 
                     function(col) {countries[which(col), ]$name}) 
# Just make sure you're correctly placing your polygon objects and point objects. 
# In this case, my point object is sp_points, and my polygon object is countries

sp_points
#> Simple feature collection with 4 features and 2 fields
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: 5413377 ymin: 4321480 xmax: 5732135 ymax: 4670058
#> CRS:            +proj=eck4 +lon_0=0 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
#>     name                geometry country_name
#> 1 Point1 POINT (5507002 4321480)  Afghanistan
#> 2 Point2 POINT (5656241 4438573)  Afghanistan
#> 3 Point3 POINT (5413377 4554775) Turkmenistan
#> 4 Point4 POINT (5732135 4670058)  Afghanistan

Old answer below, which is a way to visualize the same operation, with a completely reproducible example.

library(sf)
library(tmap)
library(dplyr)
data(World)
countries = World %>%
    filter(name == "Afghanistan" | 
             name == "Turkmenistan" |
                name == "Iran") #Picking three countries in the world

pnts = data.frame(latitude = c(34,35,36,37),longitude=c(63,65,62.5,66.5)) #Picking random points

sp_points = st_as_sf(pnts,coords = c('longitude',"latitude"))#make points spatial
st_crs(sp_points)= 4326 # Give the points a coordinate reference system (CRS)
sp_points=st_transform(sp_points,crs = st_crs(countries)) # Match the point and polygon CRS

tm_shape(countries)+
  tm_borders() +
tm_shape(sp_points)+
  tm_dots(col='red',size=1.2)

#Notice that Iran does not contain a dot

newcountries=st_intersects(countries,sp_points) #preforms the intersection
newc_logical=lengths(newcountries)>0 #Finds only countries where intersection is TRUE
newcountries2=countries[newc_logical,] #Applies it to the original countries
#The three above lines were found from: 
#https://cran.r-project.org/web/packages/sf/vignettes/sf3.html


tm_shape(newcountries2)+
  tm_borders() #Leaves you with the df newcountries2 without the country containing a dot

To answer your question about 4326, it is the coordinate reference system. This coordinate reference system is probably the most standard one for mapping, most GIS generally revert to it. The number is specifically called an EPSG code.

What you should focus on is getting both the polygon and points you are working with on the same coordinate reference system, and make sure that they are both spatial features (hence sf). The code I have showed is the most simple version of an intersection. Hopefully your problem will be simple enough to be handled this way.