I am trying to solve 3 simultaneous ODE's. Which method to use?

I am trying to solve 3 simultaneous ODE's. Which method should be used to solve these equations?? I have tried to solve it using ode45 method but I'm unable to get the answer.

dx1/df = F( f, x1, x2, x3 ) ]
dx2/df = F( f, x1, x2, x3 ) 
dx3/df = F( f, x1, x2, x3 ) 
f = [0.5 1] 
x1(0.5)=0.185,  x2(0.5)= 0.285, x3(0.5)= 0.53 
gamma = 0.577; 
xo1= 0.185; xo2= 0.285; xo3 =0.53;
theta=0.5; 
q1=1; q2 = 0.317; q3 = 0.065; 
F(x1,f) = ((q1)*(x1-gamma*(x1*f-(xo1*(1-theta)))/(f-(1-theta)))-x1*(x1-gamma*(x1*f-(xo1*(1-theta)))/(f-(1-theta))+q2*(x2-gamma*(x2*f-(xo2*(1-theta)))/(f-(1-theta)))+q3*(x3-gamma*(x3*f-(xo3*(1-theta)))/(f-(1-theta))))*f)/(((x1-gamma*(x1*f-(xo1*(1-theta)))/(f-(1-theta)))+q2*(x2-gamma*(x2*f-(xo2*(1-theta)))/(f-(1-theta)))+q3*(x3-gamma*(x3*f-(xo3*(1-theta)))/(f-(1-theta))))*f); 
F(x2,f) = ((q2)*(x2-gamma*(x2*f-(xo2*(1-theta)))/(f-(1-theta)))-x2*(x1-gamma*(x1*f-(xo1*(1-theta)))/(f-(1-theta))+q2*(x2-gamma*(x2*f-(xo2*(1-theta)))/(f-(1-theta)))+q3*(x3-gamma*(x3*f-(xo3*(1-theta)))/(f-(1-theta))))*f)/(((x1-gamma*(x1*f-(xo1*(1-theta)))/(f-(1-theta)))+q2*(x2-gamma*(x2*f-(xo2*(1-theta)))/(f-(1-theta)))+q3*(x3-gamma*(x3*f-(xo3*(1-theta)))/(f-(1-theta))))*f); 
F(x3,f)= ((q3)*(x3-gamma*(x3*f-(xo3*(1-theta)))/(f-(1-theta)))-x3*(x1-gamma*(x1*f-(xo1*(1-theta)))/(f-(1-theta))+q2*(x2-gamma*(x2*f-(xo2*(1-theta)))/(f-(1-theta)))+q3*(x3-gamma*(x3*f-(xo3*(1-theta)))/(f-(1-theta))))*f)/(((x1-gamma*(x1*f-(xo1*(1-theta)))/(f-(1-theta)))+q2*(x2-gamma*(x2*f-(xo2*(1-theta)))/(f-(1-theta)))+q3*(x3-gamma*(x3*f-(xo3*(1-theta)))/(f-(1-theta))))*f);

[Paul]

Can you show the actual code you tried? Can you at least show the definiton of the function F(f,x1,x2,x3) ?

[Dhiraj Deshmukh]

Hello Paul, I have updated the question for your reference.. kindly guide me if you have any leads.. thanks..

[Dhiraj Deshmukh]

code i have used :-

function [dy] = countermulticompfcn(f,y)
theta = 0.5;
xo1 = 0.185;
xo2 = 0.285;
xo3 = 0.53;
gamma = 0.577;
q1=1;
q2 = 0.3543;
q3 = 0.3249;
dy = zeros(4,1)
dy(1) = ((q1)*(y(1)-gamma*(y(1)*f-(xo1*(1-theta)))/(f-(1-theta)))-y(1)*(y(1)-gamma*(y(1)*f-(xo1*(1-theta)))/(f-(1-theta))+q2*(y(2)-gamma*(y(2)*f-(xo2*(1-theta)))/(f-(1-theta)))+q3*(y(3)-gamma*(y(3)*f-(xo3*(1-theta)))/(f-(1-theta)))))/(((y(1)-gamma*(y(1)*f-(xo1*(1-theta)))/(f-(1-theta)))+q2*(y(2)-gamma*(y(2)*f-(xo2*(1-theta)))/(f-(1-theta)))+q3*(y(3)-gamma*(y(3)*f-(xo3*(1-theta)))/(f-(1-theta))))*f)
dy(2) = ((q2)*(y(2)-gamma*(y(2)*f-(xo2*(1-theta)))/(f-(1-theta)))-y(2)*(y(1)-gamma*(y(1)*f-(xo1*(1-theta)))/(f-(1-theta))+q2*(y(2)-gamma*(y(2)*f-(xo2*(1-theta)))/(f-(1-theta)))+q3*(y(3)-gamma*(y(3)*f-(xo3*(1-theta)))/(f-(1-theta)))))/(((y(1)-gamma*(y(1)*f-(xo1*(1-theta)))/(f-(1-theta)))+q2*(y(2)-gamma*(y(2)*f-(xo2*(1-theta)))/(f-(1-theta)))+q3*(y(3)-gamma*(y(3)*f-(xo3*(1-theta)))/(f-(1-theta))))*f)
dy(3) = ((q3)*(y(3)-gamma*(y(3)*f-(xo3*(1-theta)))/(f-(1-theta)))-y(3)*(y(1)-gamma*(y(1)*f-(xo1*(1-theta)))/(f-(1-theta))+q2*(y(2)-gamma*(y(2)*f-(xo2*(1-theta)))/(f-(1-theta)))+q3*(y(3)-gamma*(y(3)*f-(xo3*(1-theta)))/(f-(1-theta)))))/(((y(1)-gamma*(y(1)*f-(xo1*(1-theta)))/(f-(1-theta)))+q2*(y(2)-gamma*(y(2)*f-(xo2*(1-theta)))/(f-(1-theta)))+q3*(y(3)-gamma*(y(3)*f-(xo3*(1-theta)))/(f-(1-theta))))*f)
end
solver:-
clear all
for index= 1:4
    myprompt2=['Enter the initial value of Y:(',num2str(index),'):'];
    initialval(index) = input(myprompt2);
end
lower = input('Enter Lower limit of Integration ');
upper = input('Enter Upper limit of Integration ');
[F,y] = ode45(@countermulticompfcn,[lower upper],initialval)

plot(F,y)
hold on

Answer 1

[Paul]

This line:

dy = zeros(4,1);

Needs to be replaced with:

dy = zeros(3,1);

After that, note that the equation for dy(1) evaluates to NaN at f = 0.5 and to -Inf at f = 0 (I only explored between 0 and 1). In the original question it was stated that the the tspan input to ode45 is [0.5 1]. So that will cause a problem because the initial value of dy is NaN at f = 0.5 and that kills the entire solution. I ran from [0.51 1] and got a non-NaN result. In the follow up comment the user enters the tspan limits. If those limits cause the solver to hit either f = 0 or f = 0.5 exactly there will be a problem. So I guess the real question is if it's expected that the equations for dy have these singulariites and how they should be handled if the solver has to integrate through them.

[Dhiraj Deshmukh]

Thank you very much Paul. Problem is solved. Your suggestions really hepled.

Answer 2

[Steven Lord]

This documentation page includes two examples, the second of which you can use as a model for your code to solve a system of differential equations.

[Dhiraj Deshmukh]

Thank you very much for guiding.