I have created a code to solve three coupled ODE but unable to plot its nature on the graph as one of the curve is imaginary .

clc
clear all

alpha=3.55*10^-3;
sigma=3.7292*10^11;
eta=rand;
TB=3.59*10^8;
delta_t=50*10^-12;
n=1.45;
c=3*10^8;
Q=141.3227*10^-4;

syms t p(t) AL(t) AS(t) 

% --- move to here ---

real_part=randn();
imag_part=randn();
complex=exp(-(real_part+1i*imag_part))^2;
f=sqrt((n*Q)/((delta_t)^2)*c)*complex;    % <-- this parameter must be defined before ODE
r=sqrt((n*Q)/(c*TB))*complex;
% --- move to here ---


Eq1=diff(AL)==1i*sigma*p(t)*AS(t);
Eq2=diff(AS)==1i*sigma*conj(p(t))*AL(t);
Eq3=(alpha/TB)*diff(p,2)+(alpha-1i)*diff(p)-((1i*TB)/2)*p(t)==eta*AL(t)*conj(AS(t))-1i*f;
[VF,subs]=odeToVectorField(Eq1,Eq2,Eq3);
 disp(VF)

disp(subs)

% ftotal=matlabFunction(VF);


ftotal  = matlabFunction(VF, 'vars', {'t', 'Y'});
tspan   = [0 1].*10^-6;
ic      = [0 0 r 0];        % <-- 4 states requires 4 initial values (check the order)
[t,Y]    = ode45(ftotal, tspan, ic);



 plot(t, Y(:, 1), 'LineWidth', 2, 'DisplayName', 'AS');
 hold on;
plot(t, Y(:, 2), 'LineWidth', 2, 'DisplayName', 'AL');
plot(t, Y(:, 3), 'LineWidth', 2, 'DisplayName', 'p');

Warning: Imaginary parts of complex X and/or Y arguments ignored.

 xlabel('Time');
ylabel('Solution');
title('Solution of the Coupled ODE System');
 grid on;

The plot is not showing its display name and one of the curve is imaginary. I dont know how to plot all three at the same plane.

Answer 1

[Sam Chak]

Hi @Yogesh

You can plot the real part and the imaginary part separately to verify if this meets your requirements.

alpha=3.55*10^-3;
sigma=3.7292*10^11;
eta=rand;
TB=3.59*10^8;
delta_t=50*10^-12;
n=1.45;
c=3*10^8;
Q=141.3227*10^-4;

syms t p(t) AL(t) AS(t) 

% --- move to here ---

real_part=randn();
imag_part=randn();
complex=exp(-(real_part+1i*imag_part))^2;
f=sqrt((n*Q)/((delta_t)^2)*c)*complex;    % <-- this parameter must be defined before ODE
r=sqrt((n*Q)/(c*TB))*complex;
% --- move to here ---


Eq1=diff(AL)==1i*sigma*p(t)*AS(t);
Eq2=diff(AS)==1i*sigma*conj(p(t))*AL(t);
Eq3=(alpha/TB)*diff(p,2)+(alpha-1i)*diff(p)-((1i*TB)/2)*p(t)==eta*AL(t)*conj(AS(t))-1i*f;
[VF,subs]=odeToVectorField(Eq1,Eq2,Eq3);
disp(VF)

disp(subs)

ftotal  = matlabFunction(VF, 'vars', {'t', 'Y'});
tspan   = [0 1].*10^-6;
ic      = [0 0 r 0];        % <-- 4 states requires 4 initial values (check the order)
[t,Y]    = ode45(ftotal, tspan, ic);

subplot(221)
plot(t, Y(:, 1), 'LineWidth', 2, 'DisplayName', 'AS'), grid on, title('AS')
subplot(223)
plot(t, Y(:, 2), 'LineWidth', 2, 'DisplayName', 'AL'), grid on, title('AL')
subplot(222)
plot(t, real(Y(:, 3)), 'LineWidth', 2, 'DisplayName', 'p'), grid on, title('real(p)')
subplot(224)
plot(t, imag(Y(:, 3)), 'LineWidth', 2, 'DisplayName', 'p'), grid on, title('imag(p)')

[Yogesh]

Thank you @Sam Chak , you are a life saver...