Solving coupled ode in MATLAB

I am trying solve a system of coupled differential equantions (one of which is second order). I am searching for a numerical answer. However, code keeps giving me an error and I can't track it down. Would be greatful for any help! The code is below:

%function definitions
syms r(T)  p(T);
ode1 = diff(r,T,2) == E^2 / ((m^2)*(c^2)) - (1- rs/r)*(c^2 + (h^2)/(r^2));
ode2 = diff(p,T) == (1/r^2) * L;

%function setup
f1= odeToVectorField(ode1);
F1= matlabFunction(f1);
f2= odeToVectorField(ode2);
F2= matlabFunction(f2);
odes = {F1,F2};

%initial conditions
r_0= 5*rs;
r_prime_0 = 5*rs;
r_double_prime_0 = 5*rs;
p_0 = 0;
p_prime_0 = PI/100;
conditions = [r_0;r_prime_0;r_double_prime_0;p_0;p_prime_0];
span = [0 100];

%solver
[rsol,psol] = ode45(odes,span,conditions);

Answer 1

[Star Strider]

This should get you started:

% %function definitions
syms r(T) p(T) c E h L m rs T Y 
ode1 = diff(r,T,2) == E^2 / ((m^2)*(c^2)) - (1- rs/r)*(c^2 + (h^2)/(r^2));
ode2 = diff(p,T) == (1/r^2) * L;
% %function setup
[f1, Subs]= odeToVectorField([ode1 ode2]);
F1= matlabFunction(f1, 'Vars',{T,Y, c, E, h, L, m, rs})
% %initial conditions
r_0= 5*rs;
r_prime_0 = 5*rs;
r_double_prime_0 = 5*rs;
p_0 = 0;
PI = rand;                                              % Provide Value
p_prime_0 = PI/100;
c = rand;                                               % Provide Value
E = rand;                                               % Provide Value
h = rand;                                               % Provide Value
L = rand;                                               % Provide Value
m = rand;                                               % Provide Value
rs = rand;                                              % Provide Value
conditions = [rand(3,1)];                               % Three Equations = Three Initial Conditions
span = [0 100];
% %solver
[rsol,psol] = ode45(@(T,Y)F1(T,Y,c,E,h,L,m,rs),span,conditions);

figure
plot(rsol, psol)
grid

Note that if you want a numerical solution, the required parameter values must be defined first, before calling any of the numeric ODE solvers (such as ode45).

[Lada Nuzhna]

Thank you so so much!

[Star Strider]

As always, my pleasure!

Answer 2

[Pitt Lucht]

Hello,

I have a similar problem of three coupled differential equations. I have already adapted my code to fit the suggested code from Star Strider, but I'm still runnning into problems. The current error code I'm getting is "too many input arguments", my guess is there is a problem with the way I set up my initial conditions? If somebody has the time to have a quick look over my code it would be highly appreciated!

global kr kf lf dr dm p1 p2
 
% Parameter
kr = 1;     % formation R-centers
kf = 0.1;   % formation M-centers
lf = 100;  % loss of F-centers
dr = 0.1;   % decay R-center
dm = 1;     % decay M-center
p1 = 10000;
p2 = 0;

% function definiton
syms R(T) M(T) F(T) T Y;
ode1 = diff(R,T) == -(dr*R) + (kr*M*F);
ode2 = diff(M,T) == (dr*R) - (dm*M) + (kf*(F^2)) - (kr*M*F);
ode3 = diff(F,T,2) == (dr*R) + (2*dm*M) - (kr*M*F) - (2*kf*(F^2)) - (lf*F)+p1;

% function setup
[f1, Subs] = odeToVectorField([ode1 ode2 ode3]);
F1 = matlabFunction(f1, 'Vars',{T, Y});

% initial states
R0 = 84.99;
M0 = 1.675;
F0 = 9.975;
conditions = [R0; M0; F0];
span = [0 10];

% solver
tic;
[Rsol, Msol, Fsol] = ode23s(@(T,Y)F1(T,Y,kr,kf,lf,dr,dm,p1),span,conditions);
toc;

%Log. plotting
loglog(Rsol, Msol, Fsol);
grid on;

[Star Strider]

This should actually be a new Question.

That aside, this:

% initial states
R0 = 84.99;
M0 = 1.675;
F0 = 9.975;
conditions = [R0; M0; F0; 42];                                  % <— 4 DEs => 4 ICs!  Added ‘42’ To Test Code
span = [0 10];
% solver
tic;
[Rsol, Msol] = ode23s(F1,span,conditions);                      % All Arguments Already Included In ‘F1’ After ‘subs’ Call
toc;
%Log. plotting
figure                                                          % Using 'log' Scale Not Appropriate With Negative Values, Use 'yyaxis' Instead
yyaxis left
plot(Rsol, Msol(:,[1 3 4]));                                    % 2 Arguments Only!
yyaxis right
plot(Rsol, Msol(:,2))
grid on;
legend(string(Subs), 'Location','NW')

works.

[Pitt Lucht]

Thank yo so much! And sorry for not opening a new question.

[Star Strider]

My pleasure!

A Vote for my Answer would be appreciated!

[Pitt Lucht]

It seems to not be possible to vote for comments, therefore I voted for your original answer to Lada Nuzhna since that helped out a lot too

[Star Strider]

Thank you!

Yes, voting for an Answer is the only option in this instance.