# I have solved a set of three coupled odes in MATLAB. I have to plot another function which is a function of one of these odes. How this can be done?

4 views (last 30 days)

Show older comments

I was solving the coupled rate equations of a laser by using ODE for obtaining the fluctuations of carrier density, photon number and phase. Now I have to plot a new function for frequency fluctuation which is a function of the ode for phase. Should I include the new function in the function used for defining the set of ODES or just ouside the function?.

The function defined for coupled rate equations of the laser is:

function dy = rate_eq(t,y)

%ode of carrier density

dy(1) = (I/q) - y(1)/te - (g*(y(1)-n0)/(1+(eps*y(2))))*y(2)+((((2*y(1)/te)*(1/delt))^0.5)*1.8339)-((((2*Betasp*y(1)*y(2))/te)*(1/delt)^0.5)*0.5377);

%ode for photon number

dy(2) = (g*(y(1)-n0)/(1+(eps*y(2))))*y(2)- y(2)/tp + Betasp*y(1)/te+((((2*Betasp*y(1)*y(2)/te)*(1/delt))^0.5)*0.5377) ;

%ode for phase

dy(3) = ((a/2)*g)*(y(1)-n_bar)+ (((Betasp*y(1)/2*te*y(2))*(1/delt))^0.5)*(-2.2588);

end

%The frequency fluctuation(deln(t)) is related to the phase ode as

deln=(1/2*pi)*dy(3)

I have to plot deln vs time graph. Kindly help on this

##### 1 Comment

Askic V
on 3 Mar 2023

Please have a lookt at this example:

https://www.mathworks.com/matlabcentral/answers/638120-help-using-ode45-to-solve-3rd-order-ode

### Accepted Answer

Davide Masiello
on 3 Mar 2023

### More Answers (1)

Sam Chak
on 4 Mar 2023

Hi @Rakhi

Two methods are shown in the example. The 1st method is by logical intuition. If deln depends dy(3), and the equation for dy(3) is known, then it's a direct Plug-and-Plot approach.

The 2nd method employs the deval() function that evaluates differential equation solution structure and extracts the 1st derivative of the numeric solution produced by the ode45 solver.

%% Method 1: manually inserting the equation for dy(3)

[t, y] = ode45(@rate_eq, [0 3], [1 0 3]);

figure(1), subplot(211)

plot(t, y, 'linewidth', 1.5), grid on, legend('show', 'location', 'best'), title('Method 1')

subplot(212)

a = 1;

g = 1;

n_bar = 1;

Betasp= 1;

te = 1;

delt = 1;

dy3 = ((a/2)*g)*(y(:,1) - n_bar) + (((Betasp*y(:,1)/2*te.*y(:,2))*(1/delt)).^0.5)*(-2.2588);

deln1 = (1/2*pi)*(dy3);

plot(t, deln1, 'color', '#5f86dc', 'linewidth', 1.5), grid on, legend('deln', 'location', 'best'), xlabel('t')

%% Method 2: using deval() function

sol = ode45(@rate_eq, [0 3], [1 0 3]);

tt = linspace(0, 3, 31);

[yo, yp] = deval(sol, tt); % yo is y-out, yp is y-prime (y')

figure(2), subplot(211)

plot(tt, yo, 'linewidth', 1.5), grid on, legend('show', 'location', 'best'), title('Method 2')

subplot(212)

deln2 = (1/2*pi)*yp(3,:);

plot(tt, deln2, 'color', '#5f86dc', 'linewidth', 1.5), grid on, legend('deln', 'location', 'best'), xlabel('t')

function dy = rate_eq(t, y)

I = 1;

q = 1;

te = 1;

g = 1;

n0 = 1;

delt = 1;

Betasp= 1;

tp = 1;

a = 1;

n_bar = 1;

dy = zeros(3, 1);

%ode of carrier density

dy(1) = (I/q) - y(1)/te - (g*(y(1) - n0)/(1 + (eps*y(2))))*y(2) + ((((2*y(1)/te)*(1/delt))^0.5)*1.8339) - ((((2*Betasp*y(1)*y(2))/te)*(1/delt)^0.5)*0.5377);

%ode for photon number

dy(2) = (g*(y(1) - n0)/(1 + (eps*y(2))))*y(2) - y(2)/tp + Betasp*y(1)/te + ((((2*Betasp*y(1)*y(2)/te)*(1/delt))^0.5)*0.5377);

%ode for phase

dy(3) = ((a/2)*g)*(y(1) - n_bar) + (((Betasp*y(1)/2*te*y(2))*(1/delt))^0.5)*(-2.2588);

end

### See Also

### Products

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!