Range response for Ultrasound tissue

Hi,

I'm using the MATLAB range response guide to plot a simualtion of pulses from different boundaries in human tissue. https://uk.mathworks.com/help/phased/ref/phased.rangeresponse-system-object.html

From this I'm trying to show each line is a diffrent boundary. So at 2cm is skin, 5cm is muscle and 11cm is bone. However, I can't seemed to get the image right. The code is shown below:

fc = 5e6;                       %Operating frequency
fs = 20e6;                      %Sampling frequency
propv = 1540;                   %propgataion velocity of soft tissue
propv1 = 1580;                  %propagation velocity of muscle;
propv2 = 1580;                  %propagation velocity of bone;
c = 1540;
pri = 100e-9;                   %pulse repetition interval;
prf = 1/pri; 


Numtgts = 3;
tgtpos = zeros(Numtgts);
tgtpos(1,:) = [0.02 0.05 0.11];
tgtvel = zeros(3,Numtgts);
tgtvel(1,:) = [0 0 0];
tgtrcs = db2pow(10)*[1 1 1];
tgtmotion = phased.Platform(tgtpos,tgtvel);
target = phased.RadarTarget('PropagationSpeed',c,'OperatingFrequency',fc,'MeanRCS',tgtrcs);
radarpos = [0;0;0];
radarvel = [0;0;0];
radarmotion = phased.Platform(radarpos,radarvel);
txantenna = phased.IsotropicAntennaElement;
rxantenna = clone(txantenna);

bw = fs/2;

waveform = phased.LinearFMWaveform('SampleRate',fs, ...
    'PRF',prf,'OutputFormat','Pulses','NumPulses',1,'SweepBandwidth',fs/2, ...
    'DurationSpecification','Duty cycle','DutyCycle',0.02);

sig = waveform();
Nr = length(sig);
bwrms = bandwidth(waveform)/sqrt(12);
rngrms = c/bwrms;


peakpower = 20;
txgain = 40;
transmitter = phased.Transmitter( ...
    'PeakPower',peakpower, ...
    'Gain',txgain, ...
    'InUseOutputPort',true);
radiator = phased.Radiator( ...
    'Sensor',txantenna, ...
    'PropagationSpeed',c, ...
    'OperatingFrequency',fc);

channel = phased.FreeSpace( ...
    'SampleRate',fs, ...    
    'PropagationSpeed',c, ...
    'OperatingFrequency',fc, ...
    'TwoWayPropagation',true);

collector = phased.Collector( ...
    'Sensor',rxantenna, ...
    'PropagationSpeed',c, ...
    'OperatingFrequency',fc);
rxgain = 40 ;
noisefig = 10;
receiver = phased.ReceiverPreamp( ...
    'SampleRate',fs, ...
    'Gain',rxgain, ...
    'NoiseFigure',noisefig);

Np = 128;                                       %Building the data cube
cube = zeros(Nr,Np);
for n = 1:Np
    [sensorpos,sensorvel] = radarmotion(pri);
     [tgtpos,tgtvel] = tgtmotion(pri);
    [tgtrng,tgtang] = rangeangle(tgtpos,sensorpos);
    sig = waveform();
    [txsig,txstatus] = transmitter(sig);
    txsig = radiator(txsig,tgtang);
    txsig = channel(txsig,sensorpos,tgtpos,sensorvel,tgtvel);    
    tgtsig = target(txsig);   
    rxcol = collector(tgtsig,tgtang);
    rxsig = receiver(rxcol);                    %Receive signal
    cube(:,n) = rxsig;                          %Loading data to cube        
end


imagesc([0:(Np-1)]*pri*1e6,[0:(Nr-1)]/fs*1e6,abs(cube))
xlabel('Slow Time {\mu}s')
ylabel('Fast Time {\mu}s')


matchingcoeff = getMatchedFilter(waveform);
ndop = 128;
rangeresp = phased.RangeResponse('SampleRate',fs,'PropagationSpeed',c);
[resp,rnggrid] = rangeresp(cube,matchingcoeff);
imagesc([1:Np],rnggrid,abs(resp))
xlabel('Pulse')
ylabel('Range (cm)')