Is k-means clustering supposed to behave like this? First of all, data I'm using can be found here. 
My code is:
library(readxl)
prouni <- read_excel("C:/Users/e270860661/Downloads/prouni.xlsx", 
                 sheet = "Formatados_para_rodar")
View(prouni) 
str(prouni)

coercivo <- data.frame(prouni$mensalidade, prouni$Medicina, 
prouni$nota_integral_ampla)
summary(coercivo)
coercivo$dropador <- complete.cases(coercivo)
summary(coercivo$dropador)
final <- coercivo[coercivo$dropador == TRUE,]
final$dropador <- NULL


set.seed(100)
analise <- kmeans(final, centers = 4)
str(analise)

plot(final, col = analise$cluster)

The plot I get is: 
For context, "mensalidade" means "tuition", "medicina" is a dummy variable (1 for a given program being Medicine and 0 for not being Medicine) and "nota_integral_ampla" is equivalent to required SAT score to be approved in the program.
My problem is that clustering doesn't seem to be working "multivariably". The algorithm seems to have chosen tuition thresholds and classified observations considering only these thresholds. Is my intutition right or is kmeans supposed to do this? Is there a coding error?
I'm an economist by training so this is all very new to me, sorry if it's a poorly elaborated question.
 A: I think what is seen is an artefact of having data in different scales, using somewhat inappropriate $k$ and potentially employing an inappropriate clustering algorithm.
As the data in the sample's dimensions are on substantially different scales, it is problematic to use the raw data directly. That is because simple Euclidean distances between data-points are potentially irrelevant if the differences in scale do not translate to difference in importance. If we think that all features are equally important the easiest way to achieve that is by using normalised data in the sense of : newfinal <- apply( final, 2, scale) where the values of the features contained on the columns of newfinal are set to have mean $0$ and variance $1$. That said, I do not think that using $k=4$ is a good choice, even with standardised data. I say that based on how the data look as well as what standard metrics for choosing $k$ like the GAP-statistics (available in cluster::clusGap) suggest. Running something like: clusGap( newfinal, kmeans, K.max = 6, B = 100) strongly suggests that $k=4$ is overkill and potentially $k=2$ is a much more reliable options. You can use other metrics too (e.g. fitting a Gaussian Mixture Model and examining the relevant AIC, BIC values) but I suspect these will also suggest $k<4$. 
Check out MATLAB's tutorial on Cluster Using Gaussian Mixture Models and its relevant links. It is well written and easy to follow.  $k$-means is a special case of a GMM so this can help your understanding a lot.
Aside using $k$-means I would suggest look into a density-based approach like DBSCAN. This algorithm can be easily used by using the function dbscan::dscan. The DBSAN results suggest having two major clusters and 2-3 smaller/outlier clusters. The relevant command will be something like : analise2 <- dbscan( newfinal, eps = 1). Notice DBSCAN requires to define a smallest relevant "neighbourhood size, the parameter eps, there is very good thread here explain the role of the parameter eps in more detail.
A: Welcome to Cross Validated!
Notice that your two quantitative variables are on different scales, with mensalidade ranging from zero to ten thousand and nota_integrada ranging from 450 to 800 (say).  Since the k-means minimizes a sum of squares involving the three variables, mensalidade's variability is swallowing the variability of the other two.
Standardize both variables and give it another try. (Notice that medicina will range between 0 and 1, there is no need to standardize it too).
