I wrote this code but it takes a lot time to run even though it's not running.

clc;
clear;
ti = 0;
tf = 70E-6;
tspan=[ti tf];
k = 7E-6;
h = 1E-2;
for j = 1:100
    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));
        rand(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-9;
    o = sort(10e2*rand(1,20),'ascend');
    [T,Y]= ode45(@(t,y) rate_eq(t,y,yita_mn,N,o,k),[0:1E-8:70E-6]./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(j)
    M = (abs(r));
    plot(T,M)
end




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.80;
a = 0;
T = 1500;
tp = 1E-9;
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

[Torsten]

I wrote this code but it takes a lot time to run even though it's not running.

Isn't that contradictory ?

And what is your question ?

[SAHIL SAHOO]

I think the code is fine but it's not running, but if I change[0:1E-8:70E-6]./tp to [0:1E-10:7E-9]./tp then it's working fine, I don't know how to correct this thing.

[Torsten]

Then most probably there is a singularity in between 7E-9/tp and 70E-6/tp.

While running the code, you should print the variable t in rate_eq and see where the code gets stuck.

[SAHIL SAHOO]

can you please tell me how can i do that and if do that then how can I fix that problem and we can't run the program so it's not possible to check t, am i right? if not please help me in this to correct this.

[Torsten]

Add the line

disp(t)

right after the function line

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

When the code gets stuck at a certain time, stop the code, change tspan with a value for tend some instant before the time the code got stuck and inspect the results at tend to find the reason for failure.