how can I modify my code so that whenever my Y axis reaches 2.5E-6 it get constant, means no more variation with x axis.

ti = 0;
tf = 7E-12;
tspan=[ti tf];
k = 1E-6;
h = 1E-2;
y0= [(h)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
      (h)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
      (h)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1)); 
      (h)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1)); 
      (h)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
      (h)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
      (h)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
      (h)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1)); 
      (h)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1)); 
      (h)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
      (h)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
      (h)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
      (h)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1)); 
      (h)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1)); 
      (h)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
      (h)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
      (h)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
      (h)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1)); 
      (h)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1)); 
      (h)*rand(2,1); ((-3.14).*rand(1,1) + (3.14).*rand(1,1));
      ((-3.14).*rand(20,1) + (3.14).*rand(20,1));
      zeros(1,1);
      ];
yita_mn = [
    0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1;
    1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 
    0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0;
    0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0 0;
    0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0;
    0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0 0;
    0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0 0;
    0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0;
    0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0 0;
    0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0 0;
    0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0 0;
    0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0 0;
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0 0;
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0 0;
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1 0;
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 1;
    1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0
    ];
N = 20;
tp = 1E-12;
o = sort(10e4*rand(1,20),'ascend');
[T,Y]= ode45(@(t,y) rate_eq(t,y,yita_mn,N,o,k),tspan./tp,y0);

r = ((1/20).*( exp(i.*Y(:,3)) + exp(i.*Y(:,6)) + exp(i.*Y(:,9)) + exp(i.*Y(:,12)) + exp(i.*Y(:,15)) ...
    +exp(i.*Y(:,18)) +exp(i.*Y(:,21)) +exp(i.*Y(:,24)) + exp(i.*Y(:,27)) + exp(i.*Y(:,30)) + exp(i.*Y(:,33)) ...
    + exp(i.*Y(:,36)) + exp(i.*Y(:,39)) +exp(i.*Y(:,42)) + exp(i.*Y(:,45)) + exp(i.*Y(:,48)) + exp(i.*Y(:,51)) + exp(i.*Y(:,54))+ exp(i.*Y(:,57)) + exp(i.*Y(:,60))));



figure(1)
plot(T,Y(:,81),'linewidth',1.5)
xlabel("Time(in units of t_{p})")
ylabel("coupling")
grid on
set(gca,'fontname','times New Roman','fontsize',18,'linewidth',1.8);

function dy = rate_eq(t,y,yita_mn,N,o,k)

dy = zeros(4*N+1,1);
dGdt = zeros(N,1);
dAdt = zeros(N,1);
dOdt = zeros(N,1);
P = 0.20;
a = 0;
T = 2000;
tp = 1E-12;
Gt = y(1:3:3*N-2);
At = y(2:3:3*N-1);
Ot = y(3:3:3*N-0);
for i = 1:N
    dGdt(i) = (P - Gt(i) - (1 + 2.*Gt(i)).*((At(i)))^2)./T ; 
    dAdt(i) = Gt(i)*(At(i));
    dOdt(i) = a.*Gt(i) + o(1,i).*tp;  
    for j = 1:N
        dAdt(i) = dAdt(i) + (y(81)).*yita_mn(i,j)*(At(j))*cos(Ot(j)-Ot(i));
        dOdt(i) = dOdt(i) + (y(81)).*yita_mn(i,j)*((At(j)/At(i)))*sin(Ot(j)-Ot(i));
    end
end
dy(1:3:3*N-2) = dGdt;
dy(2:3:3*N-1) = dAdt;
dy(3:3:3*N-0) = dOdt;

n1 = (1:20)';
n2 = circshift(n1,-1);
n61 = n1 +60;
n62 = circshift(n61,-1);
n80 = circshift(n61,1);

j2 =  3*(1:20)-1;
j5 = circshift(j2,-1);
j8 = circshift(j2,-2);
j59 = circshift(j2,1);

dy(n61) = (o(1,n2).' - o(1,n1).').*tp + a.*(Gt(n2) - Gt(n1)) - (y(81)).*(y(j2)./y(j5)).*sin(y(n61)) - (y(81)).*(y( j5)./y(j2)).*sin(y(n61)) + (y(81)).*(y(j8)./y(j5)).*sin(y(n62)) + (y(81)).*(y(j59)./y(j2)).*sin(y(n80));  
dy(81) = k;

end

Answer 1

[MarKf]

For a question about a value or a plot, the code you shared is needelessly complicated, you could make it easier for the people willing to help. Then again, seeing the coding practices, I'd say you are new to this, so just a warning for next time.

Basically all of that is to do:

T = 0:0.175:7; Y=funct(T); %simplified function below

Even if some parameters change in your complex code, anyway before the plot you could simply add (since Y is just monotonically increasing):

cvalue = 2.5e-6;
yi = find(Y>cvalue);
Yconstant = Y(yi(1)); 
Y(Y>2.5e-6) = Yconstant;

which gives:

plot(T,Y,'linewidth',1.5)
xlabel("Time(in units of t_{p})")
ylabel("coupling")
grid on
set(gca,'fontname','times New Roman','fontsize',18,'linewidth',1.8);
xlim([0,8]),ylim([0,8e-6]), %added for comparison

Here is the simplified function:

function y=funct(x), y=1e-6.*x; end

[SAHIL SAHOO]

here the new Y you created, although what I need to do is that, I want to change why Y(:,81) in such maner that when it reaches to 2E-6 then it won't change any further.

[MarKf]

Like mathematically, in you formula/function? Even in that case (domain/range) this can be implemented with an if y>whatever statement, so it's basically the same. Otherwise I do not understand the question, apologies.