I am trying to solve this coupled ODE from the past few days using ode45 but there's an error at ode45 function which I am unable to understand .

close all
clear all
clc
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) f Y
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);

real_part=randn();
imag_part=randn();
complex=exp(-(real_part+1i*imag_part))^2;
f=sqrt((n*Q)/(delta_t)^2*c)*complex;
r=sqrt((n*Q)/(c*TB))*complex;

ic=[0 0 r];
tspan=[1 2];

[t, sol] = ode45(@(t,Y) ftotal(t,Y), tspan, ic);

plot(t, sol(:, 1), 'LineWidth', 2, 'DisplayName', 'Eq1');
hold on;
plot(t, sol(:, 2), 'LineWidth', 2, 'DisplayName', 'Eq2');
plot(t, sol(:, 3), 'LineWidth', 2, 'DisplayName', 'Eq3');
xlabel('Time');
ylabel('Solution');
title('Solution of the Coupled ODE System');
grid on;

I have attached the question and the code which I am using it to solve the equations but there is an error at ode45 which I cant understand.

[Ayush]

could you please paste your code in the code snippet? It is cumbersome to do it from the png.

[Sulaymon Eshkabilov]

Please post your whole code. The community can simulate your code and correct it instead of re-typing the whole code.

It looks like you have got a heavy part done.

Answer 1

[Sam Chak]

Hi @Yogesh

The first issue is that the free parameter "f" must be defined before its inclusion in the ODEs. All constants and parameters should be moved before the ODE section. The second issue is that since there are 4 states, 4 initial values are required, but only 3 initial values are provided. Theoretically, the code should run after fixing these two main issues. However, the ODE solver fails to integrate from the beginning. This is likely caused by the complex-valued ODEs. Please check the equations.

%% Constant
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) f Y

% --- move to here ---
%% Parameters
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 ---

%% ODE section
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)

VF = 

subs = 

% disp(VF)
% disp(subs)
% ftotal=matlabFunction(VF);

%% Solving the ODEs
ftotal  = matlabFunction(VF, 'vars', {'t', 'Y'});

ftotal = function_handle with value:
    @(t,Y)[conj(Y(3)).*Y(2).*3.7292e+11i;Y(1).*Y(3).*3.7292e+11i;Y(4);conj(Y(1)).*Y(2).*3.771114367965276e+10+Y(3).*1.815225352112676e+19i+Y(4).*(-3.590000000000001e+8+1.011267605633803e+11i)+-1.035299568074342e+26-9.378225606249226e+25i]

tspan   = [1 2];
ic      = [0 0 r 0];        % <-- 4 states requires 4 initial values (check the order)
ySol    = ode15s(ftotal, tspan, ic);

Warning: Failure at t=1.000000e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.552714e-15) at time t.

%% Plotting section
% plot(t, sol(:, 1), 'LineWidth', 2, 'DisplayName', 'Eq1');
% hold on;
% plot(t, sol(:, 2), 'LineWidth', 2, 'DisplayName', 'Eq2');
% plot(t, sol(:, 3), 'LineWidth', 2, 'DisplayName', 'Eq3');
% xlabel('Time');
% ylabel('Solution');
% title('Solution of the Coupled ODE System');
% grid on;

[Yogesh]

Thank you so much @Sam Chak