Formulating a function on Matlab for the Shannon interpolation formula
Given your samples $x[n]$, here's a code to perform the interpolation by the sinc-functions:

F = 1000; % the sampling frequency used for drawing the continuous function
Fs = 50; % the sampling frequency
T = 1/Fs; % the sampling time

t = 0:1/F:1; % the time-samples for "continuous" plotting
ts = 0:T:1;  % the time at the sample points

% the original function. In your case, this would be integral expression
x = @(t) sin(2*pi*10*t) + cos(2*pi*11*t); 
xn = x(ts);  % the discrete sequence 

% perform interpolation from xn to x:
interpolated = 0;
for n=1:length(ts)
    % n-1 due to the 1-based indexing of matlab
    interpolated = interpolated + xn(n) * sinc((t-(n-1)*T)/T);
end

hold off;
plot(t, x(t), 'r', 'linewidth', 3); hold on;  % plot continuous
stem(ts, xn);  hold on; % plot sampled
plot(t, interpolated, 'k-o'); % plot interpolated on top of red curve




As a bonus, here's what happens when the signal is not enough bandlimited: I changed the frequency of the first sine from $sin(2\pi 10 t)$ to $sin(2\pi 30 t)$:


Note: The Shannon-Interpolation requires bandlimited signals, so in particular it does not work exactly with time-limeted signals (since these are never band-limited). So, you will see some deviations especially at the edge of the signal.