Sampling and recovery using matlab
This line is wrong:

xr(k)=xr(k) + xb(p)*sinc(fs(i)*(k-t2(p)));


The correct one is:

xr(k)=xr(k) + xb(p)*sinc( ((k-1)*Tr-(p-1)*Ts)/Ts );


Where $Ts$ is the original sampling period and $Tr$ is the reconstruction sampling period. In your code $Ts$ is given by $T_s = 1/500$ for sampling frequency of $500$ Hz. And $T_r = 10^{-4}$ set by your continuous plot period. 

Note that you will still have slight errors due to the fact that the input signal is not being perfectly bandlimited.

Here is the resulting plot for $Fs=500$ Hz, $f_0 = 200$ Hz and $T_r=10^{-4}$ seconds:


it's ok now.