FMCW Radar - Range-Doppler Plot
So the problem is that when simulating FMCW signals, there is some ambiguity in the convention on the mixing order. MATLAB's documentation makes reference to it here:

&quot;The mix operation reverses the Doppler shift embedded in x, because
of the mixing order of xref and x. The mixing order affects the sign
of the imaginary part of the output argument, y. There is no
consistent convention in the literature about the mixing order. This
function and the beat2range function use the same convention. If your
program processes the output of dechirp in other ways, take the mixing
order into account.&quot;

So you're amost there. MATLAB's functions account for this so you need to account for the change yourself when doing the manual range-Doppler map. You can do this by applying a complex exponential to each pulse that shifts the phase in the other direction so that it will produce the right Doppler shift.
For example, in the simulation loop right before you assign the dechirped signal to xr do something like the code below. I currently don't have the Phased Array Toolbox installed so I did some of the calculations manually:
for m = 1:Nsweep

    % After all the target sim stuff do your own Doppler calculations
    rangeToTarget = norm(tgt_pos - radar_pos);
    radialVel = dot(tgt_vel - radar_vel, tgt_pos - radar_pos)/rangeToTarget;

    % Negate Doppler so that closing targets produce positive frequencies
    fd = -2*radialVel/lambda;

    % De-chirp and apply phase shift for this pulse to reverse Doppler
    % to the correct value
    xr(:,m) = dechirpsig*exp(1i*2*pi*(2*fd)*(m - 1)*tm);


end

This uses the simple &quot;stop and hop&quot; model where the target moves in discrete time steps which are equal to the sweep time tm.
Addressing Range Issue
I totally to forgot to address the range part of your question. The reason this is happening is because when you use MATLAB's functions to generate the range-Doppler map, you are giving it a sampling rate $f_s = 4 \space GHz$, which in turn yields a range-bin (or range-cell) size $\delta R $ of
$$\delta R = \frac{c}{2B} = \frac{c}{2(4e9)} = 0.0375 \space m$$
If you had set RangeFFTLength = numel(waveform) when doing your own 2D FFT, you would have a HUGE number of samples in the range dimension, but it would be correct given the achieved range-bin resolution $\delta R$ and will give you the correct target ranges.
The problem now comes that when you use less range samples than what you started with by setting RangeFFTLength = 2048. This decimates the effective sampling rate, so your new range-bin size is now going to be artificially larger. If you take the ratio of the number of samples produced by using $f_s = 4 \space GHz$ and the size of the range FFT of 2048, your effective sampling rate is now 98 times worse! So your new range bin size becomes
$$\delta R_{new} = 98\frac{c}{2B} = 98\frac{c}{2(4e9)} = 3.675 \space m$$
Applying this range-bin size to the peak range index should yield the correct range in your range-Doppler map. However, note that this is still wrong. I say that because the whole point of putting that much bandwidth into a chirp is so that each range-bin in the range-Doppler map is equal to to 0.0375 m and NOT 3.675 m.
If you did it correctly, you would be able to sample at a much lower rate since your beat frequencies after mixing have lower bandwidth requirements, giving you a smaller and reasonable number of samples. Each range-bin size is still 0.0375 m, preserving the range resolution desired.