I couldn't understand why this program took a long time to run, can anyone help me to solve the problem so that the code can run fast

ti = 0;
tf = 70E-5;
tspan=[ti tf];
KC = 1;
h = 1E-2;

for j = 1:500

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))];

O = 0.0015;
p = 0.0022;
h1 = rand(1,1).*(O)+rand(1,1).*(p);
h2 = rand(1,1).*(O)+rand(1,1).*(p);
h3 = rand(1,1).*(O)+rand(1,1).*(p);
h4 = rand(1,1).*(O)+rand(1,1).*(p);
h5 = rand(1,1).*(O)+rand(1,1).*(p);
h6 = rand(1,1).*(O)+rand(1,1).*(p);
h7 = rand(1,1).*(O)+rand(1,1).*(p);
h8 = rand(1,1).*(O)+rand(1,1).*(p);
h9 = rand(1,1).*(O)+rand(1,1).*(p);
h10 = rand(1,1).*(O)+rand(1,1).*(p);
h11 = rand(1,1).*(O)+rand(1,1).*(p);
h12 = rand(1,1).*(O)+rand(1,1).*(p);
h13 = rand(1,1).*(O)+rand(1,1).*(p);
h14 = rand(1,1).*(O)+rand(1,1).*(p);
h15 = rand(1,1).*(O)+rand(1,1).*(p);
h16 = rand(1,1).*(O)+rand(1,1).*(p);
h17 = rand(1,1).*(O)+rand(1,1).*(p);
h18 = rand(1,1).*(O)+rand(1,1).*(p);
h19 = rand(1,1).*(O)+rand(1,1).*(p);
h20 = rand(1,1).*(O)+rand(1,1).*(p);


yita_mn = [
    0 h1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 h20;
    h1 0 h2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    0 h2 0 h3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    0 0 h3 0 h4 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    0 0 0 h4 0 h5 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
    0 0 0 0 h5 0 h6 0 0 0 0 0 0 0 0 0 0 0 0 0;
    0 0 0 0 0 h6 0 h7 0 0 0 0 0 0 0 0 0 0 0 0;
    0 0 0 0 0 0 h7 0 h8 0 0 0 0 0 0 0 0 0 0 0;
    0 0 0 0 0 0 0 h8 0 h9 0 0 0 0 0 0 0 0 0 0;
    0 0 0 0 0 0 0 0 h9 0 h10 0 0 0 0 0 0 0 0 0;
    0 0 0 0 0 0 0 0 0 h10 0 h11 0 0 0 0 0 0 0 0;
    0 0 0 0 0 0 0 0 0 0 h11 0 h12 0 0 0 0 0 0 0;
    0 0 0 0 0 0 0 0 0 0 0 h12 0 h13 0 0 0 0 0 0;
    0 0 0 0 0 0 0 0 0 0 0 0 h13 0 h14 0 0 0 0 0;
    0 0 0 0 0 0 0 0 0 0 0 0 0 h14 0 h15 0 0 0 0;
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 h15 0 h16 0 0 0;
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 h16 0 h17 0 0;
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 h17 0 h18 0;
    0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 h18 0 h19;
    h20 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 h19 0
    ];


N = 20;
tp = 1E-9;
o = sort(10e2*rand(1,20),'ascend');
[T,Y]= ode45(@(t,y) rate_eq(t,y,yita_mn,N,o),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))));

M(j) = max(abs(r));
end 

n = M(M<0.9)

h1 = rand(1,1).*(O)+rand(1,1).*(p);
h2 = rand(1,1).*(O)+rand(1,1).*(p);
h3 = rand(1,1).*(O)+rand(1,1).*(p);
h4 = rand(1,1).*(O)+rand(1,1).*(p);
h5 = rand(1,1).*(O)+rand(1,1).*(p);
h6 = rand(1,1).*(O)+rand(1,1).*(p);
h7 = rand(1,1).*(O)+rand(1,1).*(p);
h8 = rand(1,1).*(O)+rand(1,1).*(p);
h9 = rand(1,1).*(O)+rand(1,1).*(p);
h10 = rand(1,1).*(O)+rand(1,1).*(p);
h11 = rand(1,1).*(O)+rand(1,1).*(p);
h12 = rand(1,1).*(O)+rand(1,1).*(p);
h13 = rand(1,1).*(O)+rand(1,1).*(p);
h14 = rand(1,1).*(O)+rand(1,1).*(p);
h15 = rand(1,1).*(O)+rand(1,1).*(p);
h16 = rand(1,1).*(O)+rand(1,1).*(p);
h17 = rand(1,1).*(O)+rand(1,1).*(p);
h18 = rand(1,1).*(O)+rand(1,1).*(p);
h19 = rand(1,1).*(O)+rand(1,1).*(p);
h20 = rand(1,1).*(O)+rand(1,1).*(p);

G = (h1+h2+h3+h4+h5+h6+h7+h8+h9+h10+h11+h12+h13+h14+h15+h16+h17+h18+h19+h20)./(20)

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

dy = zeros(4*N,1);
dGdt = zeros(N,1);
dAdt = zeros(N,1);
dOdt = zeros(N,1);
P = 0.05;
a = 1;
T = 1000;
tp =1E-9;
Gt = y(1:3:3*N-2);
At = y(2:3:3*N-1);
Ot = y(3:3:3*N-0);
k = 1E-3;
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) + yita_mn(i,j)*(At(j))*cos(Ot(j)-Ot(i));
        dOdt(i) = dOdt(i) + 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);

U = sort(rand(1,20).*(0.005)+rand(1,20).*(0.003),'ascend');

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

Answer 1

[Torsten]

The reason is this line:

U = sort(rand(1,20).*(0.005)+rand(1,20).*(0.003),'ascend');

MATLAB integrators are not suited to cope with random variables in the differential equations to be solved.

These are stochastic differential equations, and the methods to solve them fundamentally differ from those for deterministic ODEs.