How to code an event that causes a shift in the resulting graph of a nonlinear dynamical system

I have a Matlab script that generates periodic oscillations in a nonlinear dynamical system. I am getting the oscillations perfectly. But when I try to insert a small disturbance (event), the oscillations are dying out. The code is shown below.

clear; clc;

Storing the rate constants/parameters using Struct function 'p'. fprintf('\nStart: %s ', datestr(datetime('now')));

start = tic;
rng('default');
NSim = 1;
[avg_p, mean_avg_vrnce] = deal(zeros(1,NSim));
[ca_, ip3_,hopen_]= deal(zeros(10001,NSim));
start1 = tic;
for A=1: NSim
    %%standard values
    p = struct('c0',2,    'c1',0.185,   'v1',6,      'v2',0.11,  'v3',0.9,...
        'k3',0.1,  'd1',0.13,    'd2',1.049,  'd3',0.943,  'd5',0.082,...
        'a2',0.2,  'tmax',100,   'ipr',10000,  'dt',0.01,  'alpha',0.2,...
        'k4',1.1,  'V4',1.2,     'Ir',1);
      %simulation time
      t = (0: p.dt : p.tmax)';
      %total time steps
      t_count = length(t) - 1;
      % Standard value
      %Initial Calcium and IP3 concentrations in microMols
      [ca] = deal([.1; zeros(t_count, 1)]);
      [ip3] = deal([0.35; zeros(t_count, 1)]);
      fprintf('\nRunning %d iteration\n', A);
      %other parameters
      alpha_h = p.a2*p.d2*(ip3(1) + p.d1)/(ip3(1) + p.d3)*p.dt;
      %randomly generate states of the channels
      ini_states = unidrnd(2,p.ipr,3) - 1;
      %initializing the open probability of the channels
      hopen = [length(nonzeros(all(ini_states,2)))/p.ipr; zeros(t_count,1)];
      isopen_temp = ini_states; %storing the states to a different matrix
      for j = 1:t_count
          %  fprintf('\nInner loop %d iteration\n', i);
          m_inf_3 = (ip3(j)/(ip3(j) + p.d1))^3/p.ipr;
          n_inf_3 = (ca(j)/(ca(j)+p.d5))^3;
          J1 = (p.v1 * m_inf_3 * n_inf_3 * hopen(j) + p.v2) * (p.c0 - (p.c1+1)*ca(j));
          J2 = p.v3*ca(j)^2/(p.k3^2 + ca(j)^2);
          ca(j+1) = ca(j) + (J1-J2)*p.dt;
          beta_h = p.a2*p.dt*ca(j); %closing rate of IPR
          ip3(j+1) = ip3(j) + (((ca(j) + p.alpha*p.k4)./(ca(j) + p.k4) ).*p.V4 - p.Ir .* ip3(j))*p.dt;
          y = rand(p.ipr,3);  %probability of changing the states
          indces = (isopen_temp & y<beta_h) | (~isopen_temp & y<alpha_h);
          isopen_temp(indces) = ~isopen_temp(indces);
          hopen(j+1) = length(nonzeros(all(isopen_temp,2)));

          if (j >= 5000 ) % disturbance at this time instant
              s = 0.05;
              ca(j+1) = ca(j) + ((ca(j+1)-ca(j))+(J1-J2))*p.dt* s;
              ip3(j+1) = ip3(j) + (((ca(j) + p.alpha*p.k4)./(ca(j) + p.k4) ).*p.V4+s - p.Ir .* ip3(j))*p.dt* s;
              y = rand(p.ipr,3);
              indces = (isopen_temp & y<beta_h) | (~isopen_temp & y<alpha_h);
              isopen_temp(indces) = ~isopen_temp(indces);
              hopen(j+1) = length(nonzeros(all(isopen_temp,2)));
          end
      end
      ca_(:,A)= ca;
      ip3_(:,A)= ip3;
      hopen_(:,A)= hopen;
      [pks,tme]=findpeaks(ca,'MinPeakHeight',max(ca)/2); %obtain the peaks with their corresponding time instants
      avg_p(:,A) = sum(diff(tme)*0.01)/(length(pks)-1); %average period of oscillation
      %to obtain the variances
      mean_avg_vrnce(:,A) = var( diff(tme)*0.01);

end
m_avg_p = mean(avg_p);
m_mean_avg_vrnce = mean(mean_avg_vrnce);
%%Figures
figure();
plot(t,ca,'k',t(tme),pks,'or','LineWidth',1.5)
set(gca,'FontWeight','bold')
xlabel('Time (s)','FontWeight','bold');
ylabel('[B] (\muM)','FontWeight','bold');
text('str',['Avg. period = ',num2str(mean(avg_p)),' s',',     No. of oscillations = ',num2str(length(pks))],'Position',[35.2053285725084 0.0327449926039681 0]);
figure();
plot(t,ip3,'k','LineWidth',1.5)
set(gca,'FontWeight','bold')
xlabel('Time (s)','FontWeight','bold');
ylabel('[A] (\muM)','FontWeight','bold');
fprintf('\nFinish: %s \n ', datestr(datetime('now')))
toc(start)

To see the oscillations without the disturbance , comment the if loop.

I would like to get a result like this

where, the dashed line is the shifted oscillation (desired result) and the dotted line is the natural oscillation without the disturbance.

But, instead, I get this

Can anyone tell me how I should code in order to get the desired results?