I couldn't understand why this program is not working, is my all variable algorithm is correct?

clc
ti = 0;
tf = 10E-4;
tspan=[ti tf];
y0=[9;9;1;1;0].*10E-2;



I =0:0.1 :10 ;
Vf = zeros(length(I),1) ; 
for i = 1:length(I)
    [T,Y]= ode45(@(t,y) rate_eq(y,I(i)),tspan,y0);
    Vf(i) = (2.*(((mean(Y(:,1).*Y(:,2).*cos(Y(:,5)))).^2 + (mean(Y(:,1).*Y(:,2).*sin(Y(:,5)))).^2).^0.5))/(mean(Y(:,1).^2)+mean(Y(:,2).^2));
end
disp(Vf)

       NaN
    0.3065
    0.3449
    0.3443
    0.3617
    0.3454
    0.3619
    0.3656
    0.3589
    0.3651
    0.3622
    0.3700
    0.3750
    0.3650
    0.3792
    0.3500
    0.3687
    0.3655
    0.3635
    0.3725
    0.3723
    0.3652
    0.3573
    0.3717
    0.3687
    0.3656
    0.3690
    0.3614
    0.3664
    0.3588
    0.3699
    0.3821
    0.3566
    0.3635
    0.3592
    0.3654
    0.3608
    0.3605
    0.3619
    0.3630
    0.3488
    0.3517
    0.3569
    0.3543
    0.3494
    0.3595
    0.3527
    0.3370
    0.3590
    0.3579
    0.3559
    0.3545
    0.3431
    0.3577
    0.3573
    0.3474
    0.3441
    0.3554
    0.3269
    0.3347
    0.3372
    0.3471
    0.3389
    0.3368
    0.3298
    0.3292
    0.3219
    0.3263
    0.3443
    0.3386
    0.3260
    0.3428
    0.3151
    0.3296
    0.3156
    0.3189
    0.3200
    0.3148
    0.3134
    0.3150
    0.3299
    0.3187
    0.3071
    0.3029
    0.3232
    0.3056
    0.2895
    0.2976
    0.3067
    0.2846
    0.3003
    0.2848
    0.2940
    0.2972
    0.2976
    0.3055
    0.3036
    0.2962
    0.2932
    0.2946
    0.3023

function dy = rate_eq(y,I)

dy = zeros(5,1);
o = 2E5;  %detuning frequency
tc = 30E-9;  %photon cavity lifetime
tf = 230E-6;  %fluorescence lifetime
a1 = 0.1;  %roundtrip loass
a2 = 0.1;   
P1 = 0.2;  %pumpstrength
P2 = 0.2;  
kc = 3E-3; %critical coupling strength
s = 0.17;
k = s.*kc;

%amplitude equation
dy(1) = ((y(3) - a1) * y(1) + k * y(2) * cos(y(5))) ./ tc;
dy(2) = ((y(4) - a2) * y(2) + k * y(1) * cos(y(5))) ./ tc;

%gain medium strength
dy(3) = (P1 - y(3) * (((abs(y(1)))^2 .*I) + 1)) / tf;
dy(4) = (P2 - y(4) * (((abs(y(2)))^2 .*I) + 1)) / tf;

%phase differential equation
dy(5) = o - (k / tc) * (((y(1) ./ y(2)) + (y(2) ./ y(1))) * sin(y(5)));

end

[KSSV]

Whjat do you mean by not working? It throws any error? Be specific.

[SAHIL SAHOO]

I edited the error please check now.

[SAHIL SAHOO]

never mind i got the solution

thanks

Answer 1

[KSSV]

You need to interchange the variables to the function rate_eq. Try this:

ti = 0;
tf = 10E-4;
tspan=[ti tf];
y0=[9;9;1;1;0].*10E-2;
[T,Y]= ode45(@(t,y) rate_eq(t,y),tspan,y0);
I =0:0.1 :10 ;
Vf = zeros(length(I),1) ;
for i = 1:length(I)
    [T,Y]= ode45(@(t,y) rate_eq(I(i),y),tspan,y0);
    Vf(i) = (2.*(((mean(Y(:,1).*Y(:,2).*cos(Y(:,5)))).^2 + (mean(Y(:,1).*Y(:,2).*sin(Y(:,5)))).^2).^0.5))/(mean(Y(:,1).^2)+mean(Y(:,2).^2));
end
disp(I)
function dy = rate_eq(I,y)
dy = zeros(5,1);
o = 2E5;  %detuning frequency
tc = 30E-9;  %photon cavity lifetime
tf = 230E-6;  %fluorescence lifetime
a1 = 0.1;  %roundtrip loass
a2 = 0.1;
P1 = 0.2;  %pumpstrength
P2 = 0.2;
kc = 3E-3; %critical coupling strength
s = 0.17;
k = s.*kc;
%amplitude equation
dy(1) = ((y(3) - a1) * y(1) + k * y(2) * cos(y(5))) ./ tc;
dy(2) = ((y(4) - a2) * y(2) + k * y(1) * cos(y(5))) ./ tc;
%gain medium strength
dy(3) = (P1 - y(3) * (((abs(y(1)))^2 .*I) + 1)) / tf;
dy(4) = (P2 - y(4) * (((abs(y(2)))^2 .*I) + 1)) / tf;
%phase differential equation
dy(5) = o - (k / tc) * (((y(1) ./ y(2)) + (y(2) ./ y(1))) * sin(y(5)));
end

[SAHIL SAHOO]

would you explain me, if I change my dy =zeros(5,1) to dy = zeros(4,1) still it's working, why so?

I have 5 coupled differential equation this shouldn't happen right?

[KSSV]

Though you define,dy = zeros(4,1)...you are defining

dy(5) = o - (k / tc) * (((y(1) ./ y(2)) + (y(2) ./ y(1))) * sin(y(5)));

so your dy is 5x1.