ZOH DAC amplitude formula testing
Here is the entire simulation code. 
First, to verify the spectrum response of ZOH, I replaced the DAC_output with an impulse signal. So I can directly check if the simulation confirms the theory or not. ( As ZOH is a kind of LTI filter, the spectrum of impulse response represents the frequency response of this system. You can try other input signals like the one you originally provided.)

Input test signal:


Output Spectrum: Note that theory spectrum entirely cover the simulated spectrum.    



Ps. DAC is simulated via an up sampling and a convolution with hold filter, h = ones(1, sample_hold)

A very high pseudo analog sample frequency should be given because there is no real 'analog' in computer simulation.

You can alter different Tp or Ts (not lower than analog pseudo freq.). These codes should work.

There are some rooms for improvement about the ADC process you provided in question. You can refer the my blog.

clear all;
close all;

end_time = 2e-2;
x=0.0001:0.0001:end_time; %10KHz ADC sampling 

ADC_response = sinc(x); %% Please check the frequency response!!

fb=1e3;  %1Khz frequency

f_noise=4e3;  %1Khz frequency

yb=sin(2*pi*fb*x); %1Khz frequency sine input signal

yf_noise=sin(2*pi*4*fb*x); %4Khz frequency sine input noise

y_in=yb+yf_noise;%total signal going into the ADC

ADC_out=conv(y_in,ADC_response);%signal after ADC

%% Started here.

% ADC_out = ADC_out(1:length(y_in)); % trim for same length

ADC_out = zeros(1,length(y_in));
ADC_out(length(y_in)/2) = 1; % impulse

%% pseudo analog frequency &amp; ZOH simulation
% parameters
analog_sampling_rate = 1e8; % 100 MHz.
digital_sampling_rate = 1e4; % according to question.

%% Important Time Parameters
Ts = 1/digital_sampling_rate; % 1e-4 as given.
Tp = 2e-5; % ZOH time in analog unit (sec), tunnable. 
hold_sample = round(Tp*analog_sampling_rate); 
%  This indicates how many samples need to be hold in this pesudo analog simulation.

%% DAC Processing 
t = 0:1/analog_sampling_rate:end_time-1/analog_sampling_rate; % pseudo analog time.

% upsampling
y_up = upsample(ADC_out,analog_sampling_rate/digital_sampling_rate);

% zero order hold filter
h = ones(1,hold_sample); % zero order hold

% filter your signal
y_zoh = filter(h,1,y_up);
y_zoh(1:floor(mean(grpdelay(h))))=[];
y_zoh = [y_zoh zeros(1,floor(mean(grpdelay(h))))];

%% Time domain plot
figure;
plot(x,ADC_out);
hold on;
plot(t,y_zoh);    
legend('yb','y_zoh');

%% Freq response of discrete time signal and freq unit conversion
d_N = length(ADC_out);
d_F = fft(ADC_out,d_N);
d_freq_step = digital_sampling_rate/d_N;
d_freq = -digital_sampling_rate/2:d_freq_step:digital_sampling_rate/2-d_freq_step;
d_freq = d_freq(d_N/2+1:end)';

%% Freq response of psudo analog time signal and freq unit conversion
c_N = length(y_zoh);
c_F = fft(y_zoh,c_N);
c_freq_step = analog_sampling_rate/c_N;
c_freq = -analog_sampling_rate/2:c_freq_step:analog_sampling_rate/2-c_freq_step;
c_freq = c_freq(c_N/2+1:end)';

%% Theory transform of DAC Transfer function.

c_thy_F = abs(sinc(c_freq*Tp)*Tp/Ts); %% For verification.

%% Plot Freq domain comparison.

figure;
plot(d_freq,abs(d_F(1:d_N/2))/d_N);
xlabel('Freq. (Hz)');

figure;
plot(c_freq,abs(c_F(1:c_N/2))*d_N/c_N);
hold on;
plot(c_freq,c_thy_F);
xlabel('Freq. (Hz)');
legend('simulated','theory');