Matlab FM Demodulation and Get Rid of Phase Folding Effect

I have a matlab code to Frequency modulation and demodulation a signal. My code is work well for modulation part. My message signal is m and modulated signal is u, code plot the message signal and its integral in one graph for plotting 1. Then signal modulated with carrier and program plots the modulated signal in time domain for plotting 2. After that, by the help of some code blocks program find the frequency spectrum of modulated signal and message signal to plot graph of them for plotting 3. In demodulation part program make some fundamental calculation for FM detection, then to obtain message signal it uses filter. Last part program plots the graph of recovered signal with message signal to compare them. I summarized all code because ı do not know whre is the problem. My problem about plotting 3 when I make zoom graph 3 I see some phase foldings or like it.Graph is not symmetric according to y-axis. I did not solve this problem, I research about them and I decided to use unwrap(). Although I tried a lot, I could not be successful. How can I get rid of this phase folding with unwrap() function. Thank you. My matlab code is ;

ts = 0.0001;% Sampling interval
t0 = 0.15;  % Duration
t = 0:ts:t0;% define time vector

%% OTHER PARAMETERS 
fc = 200;  % Carrier signal frequency 
kf =50;    % Frequency deviation constant
fs = 1/ts; % Sampling frequency 

%% MESSAGE SIGNAL SIMPLY
m = 1*(t<t0/3)-2*(t<2*t0/3).*(t>=t0/3);

%% Integration of m
int_m(1) = 0;
for k =1:length(m)-1
int_m(k+1) = int_m(k) + m(k)*ts;
end

%% PLOTTING 1
figure; subplot(211);       % Message signal
plot(t,m);grid on;xlabel('time');ylabel('Amplitude');
title('m(t)');ylim([-3 2]);

subplot(212);plot(t,int_m);% Integral of message signal
grid on; xlabel('time');ylabel('Amplitude');title('integral of m(t)');
ylim([-0.07 0.07]);

%% FM MODULATED SIGNAL
u = cos(2*pi*fc*t + 2*pi*kf*int_m); 

%% PLOTTING 2
figure; plot(t,u); % Modulated signal in time domain
grid on;xlabel('time');
ylabel('Amplitude');title('FM :u(t)');
ylim([-1.2 1.2]);

%% FINDING FREQUENCY SPECTRUM AND PLOTTING 3
% Frequency spectrum of m(t)

f=linspace(-1/(2*ts),1/(2*ts),length(t));
M=fftshift(fft(m))./length(t);         % Taking fourier transform for m(t) 
U=fftshift(fft(u))./length(t);         % Taking fourier transform for u(t)


figure;subplot(211); % Frequence spectrum of m(t)
plot(f,abs(M)); grid;
xlabel('Frequency in Hz');xlim([-500 500]);
ylabel('Amplitude');title('Double sided Magnitude spectrum of m(t)');

subplot(212);plot(f,abs(U)); % Frequency spectrum of u(t)
grid;xlabel('Frequency in Hz');xlim([-500 500]);
ylabel('Amplitude');title('Double sided Magnitude spectrum of u(t)');

%% DEMODULATION (Using Differentiator)
dem = diff(u);                 
dem = [0,dem];
rect_dem = abs(dem);

%% Filtering out High Frequencies
N = 80;        % Order of Filter
Wn = 1.e-2;    % Pass Band Edge Frequency.
a = fir1(N,Wn);% Return Numerator of Low Pass FIR filter
b = 1;         % Denominator of Low Pass FIR Filter
rec = filter(a,b,rect_dem);

%% Finding frequency Response of Signals
fl = length(t);
fl = 2^ceil(log2(fl));
f = (-fl/2:fl/2-1)/(fl*1.e-4);
mF = fftshift(fft(m,fl));              % Frequency Response of Message Signal                                      
fmF = fftshift(fft(u,fl));             % Frequency Response of FM Signal
rect_demF = fftshift(fft(rect_dem,fl));% Frequency Response of Rectified FM Signal
recF = fftshift(fft(rec,fl));          % Frequency Response of Recovered Message Signal

%% PLOTTING 4 
figure;subplot(211);plot(t,m);grid on;
xlabel('time');ylabel('Amplitude');
title('m(t)');ylim([-3 2]);

subplot(212);plot(t,rec);
title('Recovered Signal');xlabel('{\it t} (sec)');
ylabel('m(t)');grid;

My problem is in this third graph to show well I put big picture

k = -(length(X)-1)/2:1:length(X)/2;

Your k is not symmetric. If you work with symmetric k is it working fine?