How to solve differential equations for variable input parameters using for loop and ode45?

I have written code (shown below) that uses ode45 to solve 2 differential equations to determine particle number during perikinetic and orthokinetic aggregation for a given flow rate (Reynolds number) through a static mixer. The final particle count from the first equation serves as the initial condition for the second equation.

% Length

L_1 = linspace(0,Trans_L,1000);  % cm
L_2 = linspace(L_1(end,1000),L,1000);  % cm

% Time

t_1 = (L_1./v)';  % s
t_2 = (L_2./v)';  % s

% Perikinetic Section

N0_1 = [N0; 0];
[t_1,N_1] = ode45(@(t,N) PKM_DE_Initial_Fun(t,N,KA,Ortho_alpha,Vol_frac,G),t_1,N0_1);

% Orthokinetic Section

N0_2 = [0,N_1(end,1)];
[t_2,N_2] = ode45(@(t,N) PKM_DE_Initial_Fun(t,N,KA,Ortho_alpha,Vol_frac,G),t_2,N0_2);

function dNdt = PKM_DE_Initial_Fun(t,N,KA,Ortho_alpha,Vol_frac,G)

dNdt = zeros(2,1);
dNdt(1) = -KA * N(1)^2;
dNdt(2) = -(4/pi) * Ortho_alpha * Vol_frac * G * N(2);

end

Now I am trying to write code that will solve these differential equations for multiple different Reynolds numbers (Re). I have the following code to determine time for each section:

% Length

for j = 1:1:length(Re)
    if Trans_L(j) > L
       L_1(j,:) = linspace(0,L,1000); 
    else
       L_1(j,:) = linspace(0,Trans_L(j),1000); 
    end
    L2_(j,:) = linspace(Trans_L(j),L,1000);
end 

% Time

for j = 1:1:length(Re)
    t_1(:,j) = (L_1(j,:)./v(j));  % s
    t_2(:,j) = (L_2(j,:)./v(j));  % s
end

But I am unsure how to proceed from here. Should a for loop be used in the function or in the main code? How is indexing properly utilized?

Note: For the inputs to the function, G will vary with different Reynolds number values while the other inputs (KA, Ortho_alpha, Vol_frac) are constants.

Thanks!

Answer 1

[Ayush Gupta]

The code that will solve these differential equations for multiple different Reynolds numbers (Re) can be written as follows and the above script can be broken down into two functions mentionded afterwards.

for j = 1:1:length(Re)
    if Trans_L(j) > L
       L_1(j,:) = linspace(0,L,1000); 
    else
       L_1(j,:) = linspace(0,Trans_L(j),1000); 
    end
    L2_(j,:) = linspace(Trans_L(j),L,1000);
    
    t_1(:,j) = (L_1(j,:)./v(j));  % s
    t_2(:,j) = (L_2(j,:)./v(j));  % s
   
    %Call the ode solver function for each reynold number
   % CAlculate G for each specific Reynold's number since it varies with
   % Reyolds number, replace it with how it actually varies
   G = Re; % for simplication it is being equated to Re
    ode_solver_function(L_1(j,:), L2_(j,:), N0, t,N,KA,Ortho_alpha,Vol_frac,G,); 
end 
 % 
%function output = a2(L_1, L_2, v,t,N,KA,Ortho_alpha,Vol_frac,G)
function [t_1,N_1,t_2,N_2] = ode_solver_function(L_1, L_2,N0,KA,Ortho_alpha,Vol_frac,G)
   % Length

    %L_1 = linspace(0,Trans_L,1000);  % cm
   % L_2 = linspace(L_1(end,1000),L,1000);  % cm

    % Time

    t_1 = (L_1./v)';  % s
    t_2 = (L_2./v)';  % s

    % Perikinetic Section

    N0_1 = [N0; 0];
    [t_1,N_1] = ode45(@(t,N) PKM_DE_Initial_Fun(t,N,KA,Ortho_alpha,Vol_frac,G),t_1,N0_1);

    % Orthokinetic Section

    N0_2 = [0,N_1(end,1)];
    [t_2,N_2] = ode45(@(t,N) PKM_DE_Initial_Fun(t,N,KA,Ortho_alpha,Vol_frac,G),t_2,N0_2); 

end

function dNdt = PKM_DE_Initial_Fun(t,N,KA,Ortho_alpha,Vol_frac,G)
dNdt = zeros(2,1);
dNdt(1) = -KA * N(1)^2;
dNdt(2) = -(4/pi) * Ortho_alpha * Vol_frac * G * N(2);
end