I have the matlab codes to solve the stimulated brilluoin scattering equations(stokes vs z plot); but doesnt showing any plot or errors

clear all

clc

% Define the pde function

function ut = pde(t, u)

% Define problem parameters

global x1 xu n ncall gamma_k;

global coeff1 coeff2 coeff3 coeff4 TT1 TT2 omega1 omega2;

global dx1 dx2 g1 g2 v gB tau_s trise;

x1 = 0; % Lower boundary of spatial domain

xu = 15; % Upper boundary of spatial domain

n = 190; % Number of discretization points

gB = 5*10^(-11); %Brilluoin gain factor

v = 0.2e9;

tau_s = 10e-9;

z = linspace(x1,xu,n);

gamma_k = 1/tau_s ; % Some parameter

coeff1 = 2; % Coefficient 1 for outer region

coeff2 = 1; % Coefficient 2 for inner region

coeff3 = 2; % Coefficient 3 for outer region

coeff4 = 1; % Coefficient 4 for inner region

omega1 = 0; % Omega1 value for outer region

omega2 = 2*pi*12.8e9 ; % Omega2 value for inner region

dx1 = 2; % Inner region starting point

dx2 = 12; % Inner region ending point

g1 = v*gamma_k*omega1*gB; % g1 value for outer region

g2 = v*gamma_k*omega2*gB; % g2 value for inner region

TT1 = (omega2)^2 - (omega1);

TT2 = (omega2)^2 -(omega2);

trise = 0.1e-9;

Ps = 10e-3;

Pb = 5e-3;

Aeff = 50e-6;

Es = (Ps/Aeff)^(1/2);

Eb = (Pb/Aeff)^(1/2);

tstart = 0;

tend = 90e-9;

t = [tstart ,tend];

% Calculate the time-dependent shape of the pulse waveform A(t)

t1 = (t - tau_s/2) / (trise/0.45);

t2 = (t + tau_s/2) / (trise/0.45);

A_t = sqrt((tanh(t1) - tanh(t2))/2);

% Calculate Es at z = 0 using the equation

Es_z0 = (Es - Eb) * A_t + Eb;

% Define the pulse function

pulse = @(t) A_t;

% Spatial index

for ii = 1:n

% Set initial conditions

if ii == 1

y1(1) = pulse(t);

y2(1) = 5e-3;

y3(1) = 1e-11;

y4(1) = 0;

y5(1) = 0;

y6(1) = 0;

end

if ii == n

y1(n) = pulse(tend);

y2(n) = 5e-3;

end

u0 =[y1(1); y2(1); y3(1); y4(1); y5(1); y6(1)];

% Set solver options

options = odeset('RelTol', 1e-6, 'AbsTol', 1e-8); % Adjust tolerances as needed

% Initialize function call counter

ncall = 0;

[t, u] = ode23t(@pde, t, u0, options);

i1 = i(dx1<=i & i<=dx2);

i2 = i(i<dx1 | dx2<i);

coeff = zeros(n,1);

coeff(i1) = coeff2;

coeff(i2) = coeff1;

coeffsin = zeros(n,1);

coeffsin(i1) = coeff4;

coeffsin(i2) = coeff3;

omega = zeros(n,1);

omega(i1) = omega2;

omega(i2) = omega1;

g = zeros(n,1);

g(i1) = g2;

g(i2) = g1;

% Reshape u into column vectors

u = reshape(u,n,6);

y1 = u(:,1); % As

y2 = u(:,2); % Ap

y3 = u(:,3); % Aa

y4 = u(:,4); % dAa/dt

y5 = u(:,5); % Phia

y6 = u(:,6); % dPhia/dt

% BCs

y1(1) = pulse(t);

y2(n) = 5e-3;

y1x = -super(xl, xu, n, y1, 1)';

y2x = super(xl, xu, n, y2, -1)';

% Derivatives in t

y1t = y1x - y2.*y3.*coeff;

y1t(1) = 0; % boundary condition

y2t = y2x + y1.*y3.*coeff;

y2t(n) = 0; % boundary condition

y3t(i) = y4; % y3t = dAa/dt

y5t(i) = y6; % y5t = dPhia/dt y6t = d2Phia/dt^2

y4t = (-g.*y1.*y2.*coeffsin - (2*gamma_k*y4 + 2*y3.*omega.*y6) + y3.*y6.*y6 - TT.*y3);

y6t = (-g.*y1.*y2.*coeffsin + 2*omega.*(y4 + gamma_k.*y3) - 2*y6.*(y4 + y3.*gamma_k))/y3;

% Combine derivatives into ut column vector

ut = [y1t; y2t; y3t; y4t; y5t; y6t];

% Increment calls to pde

ncall = ncall + 1;

end

y1_solution = u(:,1:n)';

figure;

plot(z, y1_solution);

xlabel('Spatial Domain (z)');

ylabel('Stokes Field (y1)');

title('Stokes Field Profile');

grid minor;

end

Answer 1

[Shushant]

Hi SAMEEHA,

I understand that you are facing an issue where the function "pde" does not execute as expected when running the script. Instead of producing a graph or an error message, there is no output or indication of any error occurring. This issue arises because the function is not being called and hence is never executed.

To execute the function “pde”, kindly follow the below mentioned steps-

1. Remove the first two lines from your code that is - "clear all" and "clc".
2. Save the file as "pde.m".
3. In the command line execute

>> pde(1, 2) % modify the arguments of "t" and "u" as per your need.

Alternatively, to solve this issue add "pde(1,2)" before the definition of the function.

Refer to the following documentation for more information –

Declare function name, inputs, and outputs - MATLAB function - MathWorks India

I hope this helped in solving the issues you were facing.

Thanks,

Shushant