Matlab script: To code an event that causes a shift in the resulting graph of a nonlinear dynamical system -

i have matlab script generates periodic oscillations in nonlinear dynamical system. getting oscillations perfectly. when try insert small disturbance (event), oscillations dying out. code shown below.

clear; clc; %% % storing 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; 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 , 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 channels     ini_states = unidrnd(2,p.ipr,3) - 1;     %initializing open probability of channels     hopen = [length(nonzeros(all(ini_states,2)))/p.ipr; zeros(t_count,1)];      %%     isopen_temp = ini_states; %storing states different matrix     %%     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 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 @ 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 peaks corresponding time instants     avg_p(:,a) = sum(diff(tme)*0.01)/(length(pks)-1); %average period of oscillation     %to obtain 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 oscillations without disturbance , comment if loop.

i result

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

but, instead, this

can tell me how code 'event' in order desired results?