How do I interpret this fitted vs residuals plot? 
I don't really understand heteroscedasticity. I would like to know whether my model is appropriate or not according to this plot.
 A: As @IrishStat commented you need to check your observed values against your errors to see if there are issues with variability. I'll come back to this towards the end.
Just so you get an idea of what we mean by heteroskedasticity: when you fit a linear model on a variable $y$ you are essentially saying that you make the assumption that your $y \sim N(X\beta,\sigma^2)$ or in layman's terms that your $y$ is  expected to equate $ X\beta$ plus some errors that have variance $\sigma^2$. This practically your linear model $y = X\beta + \epsilon$, where the errors $\epsilon \sim N(0,\sigma^2)$.
OK, cool so far let's see that in code:
set.seed(1);            #set the seed for reproducability
N = 100;                #Sample size
x = runif(N)            #Independant variable
beta = 4;               #Regression coefficient
epsilon = rnorm(N);     #Error with variance 1 and mean 0
y = x * beta + epsilon  #Your generative model
lin_mod <- lm(y ~x)  #Your linear model

so right, how do my model behaves: 
x11(); par(mfrow=c(1,3));   #Make a new 1-by-3 plot
plot(residuals(lin_mod)); 
title("Simple Residual Plot - OK model")
acf(residuals(lin_mod), main = ""); 
title("Residual Autocorrelation Plot - OK model");
plot(fitted(lin_mod), residuals(lin_mod)); 
title("Residual vs Fit. value - OK model");

which should give you something like this:

which means that your residuals do not seem to have an obvious trend based on your arbitrary index (1st plot - least informative really), seem to have no real correlation between them (2nd plot - quite important and probably more important than homoskedasticity) and that fitted values do not have an obvious trend of failure, ie. your fitted values vs your residuals appear quite random. Based on this we would say that we have no problems of heteroskedasticity as our residuals appear to have the same variance everywhere. 
OK, you want heteroskedasticity though. Given the same assumptions of linearity and additivity let's define another generative model with "obvious" heteroskedasticity problems. Namely after some values our observation will be much more noisy.
epsilon_HS = epsilon;               
epsilon_HS[ x>.55  ] = epsilon_HS[x>.55 ] * 9       #Heteroskedastic errors

y2 = x * beta + epsilon_HS      #Your generative model
lin_mod2 <- lm(y2 ~x)            #Your unfortunate LM

where the simple diagnostic plots of the model:
 par(mfrow=c(1,3));   #Make a new 1-by-3 plot
 plot(residuals(lin_mod2)); 
 title("Simple Residual Plot - Fishy model")
 acf(residuals(lin_mod2), main = ""); 
 title("Residual Autocorrelation Plot - Fishy model");
 plot(fitted(lin_mod2), residuals(lin_mod2)); 
 title("Residual vs Fit. value - Fishy model");

should give something like:

Here the first plot seems a bit "odd"; it looks like we have a few residuals that cluster in small magnitudes but that is not always a problem... The second plot is OK, means we have not correlation between your residuals in different lags so we might breathe for a moment. And the third plot spills the beans: it is dead clear that as we got to higher values our residuals explode. We definitely have heteroskedasticity in this model's residuals and we need to do something about (eg. IRLS, Theil–Sen regression, etc.) 
Here the problem was really obvious but in other cases we might have missed; to reduce our chances of missing it another insightful plot was the one mentioned by IrishStat: Residuals versus Observed values, or in for our toy problem at hand: 
 par(mfrow=c(1,2))
 plot(y, residuals(lin_mod) ); 
 title( "Residual vs Obs. value - OK model")
 plot(y2, residuals(lin_mod2) ); 
 title( "Residual vs Obs. value - Fishy model")

which should give something like:

here the first plot seems "relatively OK" with only a somewhat hazy upward trend in the residuals of the model (as Scortchi mentioned see here as to why we are not worried). The second plot though exhibits this problem fully. It is crystal clear we have errors that are strongly dependent on the values of our observed values. This manifesting in issues with the coefficient of determination $R^2$ of our models at hand; eg. the "OK" model having an adjusted $R^2$ of $0.5989$ while the "fishy" one of $0.03919$. Thus we have reasons to believe model misspecification might be an issue. (Thanks to Scortchi for pointing out the misleading statement in my original answer.)
In fairness of your situation, your residuals vs. fitted values plot seems relative OK. Checking your residuals vs. your observed values would probably be helpful to make sure you are on the safe side. (I did not mention QQ-plots or anything like that as not to perplex things more but you may want to briefly check those too.)  I hope this helps with your understanding of heteroskedasticity and what you should look out for.
A: Your question seems to be about heteroscedasticity (because you mentioned it by name and added the tag), but your explicit question (e.g., in the title and) ending your post is more general, "whether my model is appropriate or not according to this plot".  There is more to determining if a model is inappropriate than assessing heteroscedasticity.  
I scraped your data using this website (ht @Alexis).  Note that the data are sorted in ascending order of fitted.  Based on the regression and upper left plot, it seems to be sufficiently faithful:  
mod = lm(residuals~fitted)
summary(mod)
# ...
# Residuals:
#   Min       1Q   Median       3Q      Max 
# -0.78374 -0.13559  0.00928  0.19525  0.48107 
# 
# Coefficients:
#   Estimate Std. Error t value Pr(>|t|)
# (Intercept)  0.06406    0.35123   0.182    0.856
# fitted      -0.01178    0.05675  -0.208    0.836
# 
# Residual standard error: 0.2349 on 53 degrees of freedom
# Multiple R-squared:  0.0008118,  Adjusted R-squared:  -0.01804 
# F-statistic: 0.04306 on 1 and 53 DF,  p-value: 0.8364


I don't see any evidence of heteroscedasticity here.  From the upper right (qq-plot), there doesn't seem to be any problems with the normality assumption either.  
On the other hand, the "S" curve in the red lowess fit (in the upper left plot), and the acf and pacf plots (at the bottom) do seem problematic.  At the far left, most of the residuals are above the gray 0 line.  As you move to the right, the bulk of the residuals drop below 0, then above, and then below again.  The result of this is that if I told you I were looking at a particular residual and that it had a negative value (but I didn't tell you which one I was looking at), you could guess with good accuracy that the residuals nearby were also negatively valued.  In other words, the residuals are not independent—knowing something about one gives you information about others.  
In addition to plots, this can be tested.  A simple approach is to use a runs test:  
library(randtests)
runs.test(residuals)
#  Runs Test
# 
# data:  residuals
# statistic = -3.2972, runs = 16, n1 = 27, n2 = 27, n = 54, p-value = 0.0009764
# alternative hypothesis: nonrandomness

The implication of this is that your model is misspecified.  Because there are two 'bends' in the relationship, you will want to add $X^2$ and $X^3$ terms to your model to account for that.  
To answer your explicit questions:  Your plot shows serial autocorrelations / non-independence of your residuals.  It means that your model is not appropriate in its current form.  
