Define formulas for ODE input

Hi all,

I would like to ask your help for a problem which is probably simple but I cannot find the proper way to do it. I have to solve a set of ODEs with the ode15s, defined as:

• d^2A/dy^2 + C*(dA/dy) - C*dB/dy = 0
• dB/dy = 1/C*f(B)

I am now however stuck on how to define these coupled ODEs so they can be solved by ode15s (or any other solver). Does anyone have a suggestion on how I can do this?

Answer 1

[Alan Stevens]

Here's one way:

% d^2A/dy^2 + C*(dA/dy) - C*dB/dy = 0
% dB/dy = 1/C*f(B)
% Turn the 2nd order ODE inot two 1st order ones
% dA/dy = v;  dv/dy = 1/C*f(B) - C*v;


    
yspan = [0 1];  % Or whatever the start and end values of y you want.

% Set initial values (arbitrary values here)
A0 = 1;
v0 = 0;
B0 = 1; 
AvB0 = [A0, v0, B0];

% Call ode15s
[y, AvB] = ode15s(@rates, yspan, AvB0);
%
% Extract results
A = AvB(:,1);
v = AvB(:,2);
B = AvB(:,3);

% Plot results
plot(y,A) % etc.

function dAvBdy = rates(~, AvB)

        A = AvB(1);
        v = AvB(2);
        B = AvB(3);
        f = sqrt(B); % Define your function of B here
        C = 1;       % Define your value of C here
        
        dAvBdy = [v;
                  1/C*f - C*v;
                  1/C*f];

end

[Jeroen]

Thanks!