chirp generation in for loop

Perhaps someone could help. I would like to generate a sinewave function with frequency modulations at given time points (timestamps, in the code). There should be the chance to set amplitude of the frequency excursions (same for every chirp, even better if they could be modified for each chirp but this might come later).

Here is what I did: I am using a function to generate the frequency vector containing frequency peaks of the desired excursion (max freq - min freq, in my example the chirp amplitude is 400Hz), then I am calculating a sinewave using the vector. I already know this does not really work because the freq gets multiplied per time. As a result the chirp amplitude increases over time (which I do not want). What I did (I do not have the code here but I can paste it later) was to insert a for loop to calculate individual sinewaves with different freq values taken from the freq vector (one every N, with N set differently depending on the "resolution" one wants to have).

The aim would be to generate a signal (sinewave) by appending each segment calculated in the for loop to a new "signal vector". I am wondering whether this could be a way to go or if there are better options (it might get very slow in case one uses high sampling rates and small segments). Also the problem could be that, without using a high "resolution" (small steps in the for loop) the assembled signal might be full of artifacts due to possible phase differences between pairs of consecutive segments.

Any suggestion on how to do this in a clean and efficient way ?

nfft=2^11;
overlap=round(0.5*nfft); 

fs = 20000; %set sampling frequency
ts=1/fs;    %calculate time unit
min=800;    %set baseline frequency (min) and peak freq excursion (max) 
max=1200;
T=10;       %max duration of sweep
t=0:ts:T;
Amp=1; % set amplitude of the signal produced

timestamps=[1 3 7]; % sec at which chirps are produced
    
    min_timevec=0;   % onset time for chirp production
    max_timevec=T;   % sec
    sigma=0.1;       % standard deviation of each gaussian peak (steepness)
    peak=10;         % set peak amplitude

f=min:(max-min)/(T/ts):max;

[spkvec,timevec,updatedpeak]=spikegauss(timestamps,fs,min_timevec,max_timevec,sigma,peak);



ax1=subplot(3,1,1);
plot(timevec,spkvec,'k'); 
ax2=subplot(3,1,2);



chirptrace=spkvec+800; %this is the frequency vector to be used to produce 
                       %chirps. It has peaks of amplitude = Amp, occurring at "spkvec" timestamps
y=Amp*sin(2*pi*(chirptrace.*t)); %attempt to generate a sinwave using those freq values as above
plot(t,chirptrace);
ax3=subplot(3,1,3);


spectrogram(y,hanning(nfft),overlap,nfft,fs,'yaxis');
ylim([0, 3])

ax=[ax1,ax2,ax3];
linkaxes(ax,'x');
% sound(y,fs);
% end



function [spkvec,timevec,updatedpeak]=spikegauss(timestamps,fs,min_timevec,max_timevec,sigma,peak)
% Generate SPKVEC time series from TIMESTAMPS. 
%
% Syntax: [spkvec,timevec,updatedpeak]=spikegauss(timestamps,srate,min_timevec,max_timevec,sigma,peak)
%
% Each spike is represented by a gaussian centered on each of the TIMESTAMPS
% SIGMA is the standard deviation
% PEAK is the value of the peak of the gaussian (use peak=0 for gaussian
% integral = 1; thus sum(spkvec) is equal to the number of spikes)

% Doubts, bugs: [e-mail]
 
timevec = min_timevec : 1/fs : max_timevec;
spike_count = histc( timestamps, timevec );
%gaussian kernel (mean=0; sigma from input)
gk_mu    = 0;      
gk_sigma = sigma;
gk_x = -10*sigma+1/fs : 1/fs : +10*sigma;
gk = 1/(sqrt(2*pi)* gk_sigma ) * exp( - (gk_x-gk_mu).^2 / (2*gk_sigma^2));
if peak==0
    gk=gk/sum(gk);
else
    gk=peak*gk/max(gk);
end
updatedpeak=max(gk);
spkvec=conv( spike_count , gk, 'same'); 
end