Generating IF signal for FMCW radar
The principle behind FMCW is that you transmit a chirped signal and receive a time delayed version of it after reflecting from a target. After mixing and filtering, the resulting signal is a sinusoid at a frequency that is a function of the target's range. This frequency is known as the "beat" frequency $f_b$. Thus, the dechirped signal will have a form of

$$x(t) = e^{j(2{\pi}f_bt + \, \phi)} = e^{j2{\pi}f_bt}\,e^{j{\phi}}$$

Where $\phi$ is a general phase term that we'll ignore for now as it won't effect the determination of the beat frequency, and therefore range. Also, let's not worry about Doppler for now...that's an additional phase term we can easily add later. We'll concentrate on the expression for the homodyned (mixed) signal.

Let the chirped signal we transmit be

$$s_{tx}(t) = e^{j\pi\frac{\beta}{\tau}t^2}$$

Where $\beta$ is the sweep bandwidth of the chirp and $\tau$ is the chirp's length, or pulse width. After reflecting from a target we receive the signal after some delay $t_d$, we have

$$s_{rx}(t) = e^{j\pi\frac{\beta}{\tau}(t - t_d)^2} = e^{j\pi\frac{\beta}{\tau}(t^2 - 2tt_d + t_d^2)}$$

After mixing $s_{rx}(t)$ with $s_{tx}(t)$, which is equivalent to a frequency shift, the higher order term containing $e^{j\pi\frac{\beta}{\tau}t^2}$ drops off and we're left with

$$x(t) = e^{j\pi\frac{\beta}{\tau}(-2tt_d + t_d^2)} = e^{-j\pi\frac{\beta}{\tau}2tt_d}\,e^{j\pi\frac{\beta}{\tau}t_d^2}$$

Now compare this with the first equation, paying attention to the first term, again ignoring the constant phase term. We can then equate the phase functions
$$-\pi\frac{\beta}{\tau}2tt_d = 2{\pi}f_bt$$

So that then we have

$$f_b = -\frac{\beta}{\tau}t_d$$

Since we know our pulse travels at the speed of light $c$, we can rewrite the target's delay in terms of range $R$ and yield the mapping between target range and its beat frequency

$$t_d = \frac{2R}{c} =&gt; f_b = -\frac{2R\beta}{c\tau}$$

So generating a dechirped signal is straight forward since it's simply a sinusoid at some beat frequency $f_b$.

Note that these equations apply only to upchirps and downchirps. The negative sign takes care of itself in either case. Triangular and more exotic chirps yield additional frequency terms but this process can be expanded to cover that as well.

To add Doppler, you can add a constant phase term that is updated as you collect pulses to form a range-Doppler map. You can actually start at zero phase for the first pulse and progress from there for the purpose of a simulation. Your additional phase term would look something like

$$e^{j2{\pi}f_dnT_c}$$

Where $n$ is the current pulse number starting at 0 and $T_c$ is your equivalent pulse repetition interval (PRI) which is akin to your sweep time for FMCW radars.

EDIT: After having some time to look at your code directly, I found a few issues.

First, you are missing a factor of two in the Doppler component of the phase.

Second, without getting into theory, your particular system supports a Doppler span wider than what you force the horizontal axis to. This will erroneously change where you perceive the target to be.

Third, the time vector you are using to generate the beat frequency must be reset to 0 after every pulse. This is because the time vector needs to be relative to the time delay of the target $t_d$.

Here is your modified code. I currently do not have the Phased Array toolbox to generate and display the range-Doppler map, so I did it manually.

c = 3e8;            %speed of light
range_max = 180;    %max detection range
tm = 6*(2*range_max/c); %sweep time
%tm is 7.2e-6 s
bw = 200e6;         %sweep bandwidth
sweep_slope = bw/tm;
v_max = 150*1000/3600;  %target max velocity
fc = 77e9;          %radar frequency
lambda = c/fc;      %radar wavelength
fs = 72e6;          %sampling rate
%sampling rate based on ADC datasheet
chirps = 64;        %frame size
samples = ceil(tm*fs);  %samples in one chirp

%% target
R0 = 20; %range in meters
V = 40; %radial velocity, m/s

%% 
t = 0; %time
mix = zeros(samples, chirps); %mixer output
for i=1:1:chirps
    td = 2 * R0 / c; %round trip delay
    phi0 = 4*pi*fc*R0/c; %inital phase
    t = 0; % Reset

    for j=1:1:samples
        a = (-2*pi*fc*2*V*i*tm/c ...    %phase shift
             -2*pi*(2*V*(fc+i*bw)/c + sweep_slope*td)*t);   %frequency
        mix(j,i) = 0.5*cos(a);
        t = t + 1/fs;
    end
end

%% Form the range-Doppler map (RDM)

% RDM axes
rangeBinAxis = (0:samples-1).*c/(2*bw);
dopplerBinSize = (1/tm)/chirps;
velocityBinAxis = (-chirps/2:chirps/2-1).*dopplerBinSize*lambda/2;

% 2D FFT to perform range and Doppler compression (i.e. form the RDM)
rdm = fftshift(fft2(mix), 2);

% Plot the RDM for the valid ranges of interest - targets ahead of you
figure;
surf(velocityBinAxis, rangeBinAxis(1:ceil(samples/2)), 20*log10(abs(rdm(1:ceil(samples/2), :))));
% surf(velocityBinAxis, rangeBinAxis, 20*log10(abs(rdm)));  % See the entire spectrum
xlabel("Range (m)");
ylabel("Velocity (m/s)");
axis tight;
shading flat;
view(0, 90);
colorbar;

% figure(1)
% rngdopresp = phased.RangeDopplerResponse('PropagationSpeed',c,...
%     'DopplerOutput','Speed','OperatingFrequency',fc,'SampleRate',fs,...
%     'RangeMethod','FFT','SweepSlope',sweep_slope,...
%     'RangeFFTLengthSource','Property','RangeFFTLength',2048,...
%     'DopplerFFTLengthSource','Property','DopplerFFTLength',256);
% 
% clf;
% plotResponse(rngdopresp,mix);
% axis([-v_max v_max 0 range_max])


Some examples to show the target being positioned appropriately in the RDM (scale is in dB):

R0 = 32 m , v = 40 m/s


R0 = 150 m , v = 40 m/s