Well, chaos and bifurcations are essentially manifested as a result of minor changes in the boundary conditions. I suggest, try changing the initial conditions and solve the equations again.

If the right-hand side of the differential equation dx/dt = f(t,x) is Lipschitz and the initial conditions are given (x(t0) = x0) then the solution is unique. Changing the initial conditions will change the solution. Numerical methods will approximate that unique solution, but if you change the initial conditions, the numerical method will approximate the solution for the new initial conditions.

4th order Runge-Kutta is a very simple integrator, for high accuracy you may need set the time step size very small, but this lead the round off error increase. Then with the same eqution, you may get different result in differnet compute for a long-term integration.

For a complex chaos system, maybe you need a better integrator, at least should have local error control, for example, 8th Dormand-prince. Then you can use bigger time step size, know  and control the local error during whole integration. And also know how much error will produce during the whole inegration.

Base on a high accurrate integrator, give some disturbance in initial condition or in some point, you can see the bifurcations. From the derivation processing of your analytic solution, your should know how the system will follow different trajectory when some disturbance introduced. And you should know how much the disturbance (or error) should be added.

I agree with Xiaolong Zhang you should change the numerical method but you should also look at implicit numerical methods such as Radau5 or BDF which are very efficient and with strong stability properties. 

As David already pointed out, for given initial conditions any proper differential equation  has a unique solution. You should tell us about the particular circumstances that let "the mathematical methods result in multiple solutions". For truely chaotic equations it is impossible to follow the 'unique exact solution' over an extended time span even if one uses the most sophisticated integrators and (an important degree of freedom not yet mentioned in the discussion) multiple precission arithmetics.

Dear researchers

I appreciate your answers. Although I agree with the presence of round off error in Runge-Kutta method, but the solution extracted is in agreement with those of the analytical method. However, the main problem is finding the other possible solutions, which are in the set of analytical results. I also tried adding disturbance to the initial conditions but the same steady state solution can be extracted.

For example in the following article, the numerical and analytical solutions of Eqs. (26) are depicted in Fig. 2. For some of the frequencies, there exist two stable and one unstable solutions.

Article Stability and bifurcation of an axially moving beam tuned to...

Dear Leon, when solving an initial condition differential equation, you do not have any initial estimation.

Dear Leon, To my knowledge there does not exist any numerical method with initial estimates to solve initial condition differential equations. If you know any, please let me know. 

@Saeed

The problem treated in Eq. 26 and Fig. 2 is not an initial value problem. Unfortunately the authors do not clearly indicate to which differential equation and to which initial conditions they applied Runge-Kutta-4 in order to obtain the red circled points in Fig. 2.

From an engineering perspective it is not a natural approach to try to understand periodically driven non-linear oscillations by studying the initial value problem. The technical mechanics literature has ingeneous methods to find out properties of periodic amplitude components which do not refer to the initial conditions which realize them. Any resonable discussion of the Duffing oscillator helps to understand the relation between the initial value problem and the properties of resonance amplitudes.

For  problems like Eq. 26 there is no need for multiple precission arithmetics or  sophisticated integrators since the appearance of multiple resonance branches has simple algebraic reasons.   

I suggest two methods:

1) try implementing adaptive RK method with varying time-steps. Have a look at the article: https://arxiv.org/pdf/1602.08652v3.pdf

2) Invoke perturbation effects as done in http://journals.aps.org/pre/abstract/10.1103/PhysRevE.56.2858

Send the differential equation in pdf format . Hope to give some idea.

Prof Rath

Dear Rath, an example of the problem is given in the answers of previous page.

@Saeed

Since you by now should know that the multiple solutions in the paper you cited do not appear in the well-defined initial value problem but in the much less defined problem of amplitude and phase responses, your answer to Professor Rath is inadequate. As it looks to me, you are not really interested in the matter. If, however, you are willing to learn what is going on in your reference paper, here is a very accessible source:

Dear Mutze, I appreciate all answers here sincerely. My problem is the same as the article I mentioned. However I have 6 differential equations in the set. Since the solution is decaying overtime, and I am looking for the steady state solution, therefore for any initial conditions I get the same steady steady solution. The function of the ODEs in Matlab is also attached. For gm = 0.0001, sigma0 = 400, f0 = 1000, there exist 2 stable and 1 unstable solution. Two steady state analytical solutions for this case are:

X1 = [0.9221 -0.00627 0.3438 -0.00619 -0.00102 0]

X2 = [-0.0881 -4.66884e-06 -8.705e-05 -7.4635e-08 5.098e-12 1.3892e-14]

@Saeed

Unfortunately, I neither see a differential equation nor a solution in your respond.

At least you give a name to the kind of solutions you are searching: 'steady state solution'. You probably are aware of the fact that there is no Runge-Kutta method for finding steady state solutions of differential equations.  Mathematica, for instance, has no algorihm for finding steady state solutions. The article concerning the Duffin equation, which was appended to my previous respond, makes the reasons for this quite clear.

If you hope to get useful information from the RG community you should invest a bit of work and bring your problem into an unambiguous form.

@Ulrich

The file FunC is the set of six first order differential equations in the state space form of DX = (AX + Nonlinear terms). The only way to find the steady state solution numerically is to find solution for 0

If you are indeed looking for steady-state solutions, then why do you need to solve a differential equation? When the derivatives are suppressed then there are six simultaneous nonlinear equations in six unknowns. These could be solved using a multidimensional Newton-Raphson method. Granted, it might be quite difficult finding an initial iterate which will converge to the third (possibly unstable in time) solution, but it will be possible if a third solution exists. When it is found, you could then check if the solution is unstable by using it (plus a very small perturbation) as an initial condition for your timestepping code.

Moreover, if the equations you are solving have parameters which you intend to vary, then such a Newton-Raphson code could be used to follow the trajectory of the steady solutions in parameter space.

Dear Rees, Thank you for your answer.

I am indeed looking for the time response, and only verifying my solutions by comparing the analytical steady-state solutions to the time response at e.g. t = 500 sec. However, I have already found some time response solutions by some perturbations in the initial conditions, but this try and error based approach is really time consuming.

@Andrew

Steady-state solution is an unsharp term here. It seems not to mean time-independent solution (i.e. stationary solution). At least the Eq. 26 to which Saeed and I referred several times has obviously no time-independent solution.

@Saeed

If you next respond fails to give an understandable description of your problem I will no longer follow this tedious discours.

@ Ulrich

Sir, I have six first order non-linear differential equations as dXn/dt = f(Xn) (for x = 1...6). In which f is a nonlinear equation and these equations are given in the FunC file. I also have the analytical steady-state solutions of these equations, e.g. for sigma0 = 400, I have 2 stable steady-state solutions (i.e. X1 and X2 in previous responses), and 1 unstable solution. Now I solve the differential equations for zero initial conditions via Runge-Kutta (as in Code file). As a result I come to 6 time-dependent solutions which are plotted when running the file Code. I check the point on these figures for e.g. t = 300, and I see that these points correspond to the X2. Now I want to solve the differential equations again and find other response that their steady-state solutions (i.e. for t =300) corresponds to X1.

I don't know how to further put my question, since the equations are long and I do not have them on PDF.

Sincierly

@Saeed

If, the right-hand side of your non-linear ODE of first order does not depend explicitely on time and, in particular, does not contain a periodic driving force term, why you described your problem to Professor Rath as analog to the Huang Su Chen paper? (This paper deals with periodically driven problems for which all solutions depend nontrivially on time).

In the framework which you now suggest the word 'response' makes no sense, since there is no stimulus.

Does 'want to solve the diffferential equations again' mean that you vary the initial conditions of your Runge-Kutta integration? If yes, according to which plan?

I understand that you refrain from reproducing the ugly numbers from your funcC.m file. But it should be possible with limited effort to make clear on which variables the right-hand side of your ODE depends. In particular: does time come in in any other form than X_i(t), i=1,...6 ? E.g. as f_0(t) ?

Given that Andrew pointed out that finding steady-state solutions is not a differential equations problem at all, and that I indicated that 'steady-state solution' may here used in a meaning different from Andrew's interpretation, I would have expected clear words from you what your definition of steady-state solution is.

There is still a bit of hope on my side that the whole discussion may make sense in the end.

@Ulrich

The right side does not have a periodic term or any time dependent term but has a constant "f0". These equations are appended in the PDF file. In this file p1, p2, p3, q1, q2, and q3 are time dependent, all alphas are constants, and alpha7r = f0. In other words f0 is the stimulus.

The problem is I do not know the plan to solve the Eqs again to find the second solution.

Andrew meant that I assume DX = 0 (dpn/dt = 0, dqn/dt = 0) and solve the resulting equations by  Newton-Raphson method. It is possible, however I want to find the whole response and only compare the steady-state solutions.

By steady state solution, I mean p1, p2, ..., q3 for ( say t = 300 sec), this steady-state solution can also be obtained by Andrew's suggestion.

At last I apologize to refer to the mentioned publication, since the equations are a bit different.

@Saeed

What's about the sigmas and the omegas?

When you write that you 'want to find the whole response' what does this mean for the given equations? As far as I'm informed the concept of response for non-linear systems tries to be analog to the corresponding well-defined concept for linear systems and is finally based on linearization of the system.

Mathematica has a package on 'control systems' which also treats non-linear system. It does so in a language of system theory with which I'm no longer really familiar. Probably MatLab has something similiar. If you would describe your concept of response I would try to figure out whether Mathematica has methods to compute that response. 

@Ulrich

sigmas are also constants, sigma0 = 400, and other ones are given in FunC file.

About the whole response, I mean not just the steady-state solution but the time-dependent functions p1(t), p2(t), ... , q3(t).

I appreciate your time and apologize for the inconvenience.

@Saeed

This can't be the understanding of 'response'. Since for any 8 numbers c1,...c6, t0, f0 you get a solution t |--> (p1(t),...,q3(t)) which satisfies the ODE and p1(t0) = c1, ... q3(t0) = c6. You need an idea on which of these many, many solutions should represent the response to the stimulus f0.

Perhaps your application is such that it makes sense to work with the initial conditions directly (e.g. since they can be set in the experiment or apparatus under consideration). Then one would not have to struggle with shaky notions ('what, after all, is a response ?') but would simply get the solution by Runge-Kutta or any other numerical solution of the initial value problem.

As it looks to me, your attempt to communicate your problem as a problem dealing with a specific ODE did not work. One would need some of its context to understand it.

Dear Saeed,

There are too many replies already for me to look back over them all, but I don't recall seeing what the application is for your six equations. However, all your right hand side terms are either linear or cubic in the dependent variables. The names of your coefficients also lead me to think that each p is the real part of a complex dependent variable and q the imaginary part. This could, of course then mean that you have three complex Landau equations which you have decomposed into their real and imaginary parts - and the equations themselves generally arise from a weakly nonlinear theory of some sort. Such equations do have a variety of potential solutions which include evolution towards a steady state (i.e. constant in time), and a finite-time blow-up. Clearly you don't have the latter, but the former appears to be the case for you.

My understanding and experience of coupled Landau equations doesn't extend to those with complex coefficients, but solutions of those with real coefficients does suggest that you will have an exponentially decaying approach to a steady state. It also suggests that there is at least one unstable steady-state solution which you won't find using an unsteady solver, but which could be found using the NR method discussed earlier. If your system has more than one steady-state (t->infinity) solution, then each will have its own domain of attraction, i.e. sets of initial conditions which will tend towards one or other of the stable steady states.

In all of this I agree with Ulrich that insufficient information was given at the outset of this discussion, and it has taken too long to get to what is important for the problem you are solving.

Not knowing exactly the problem at hand, it is hard to say. For some classes of differential equations and certain types of bifurcation, numerical continuation methods might give you some control over which solutions to track.

A problem that does not have unique solution isn't well posed. it means you do not have enough boundary conditions. you can enforce it by imposing some smoothness requirements.

Probably, by mathematical method you man analytic approaches.

I suggest First integral method for different solutions.

Maybe you may have issue with boundary conditions, that is the reason for non-uniqueness of the solution. and make sure you set you differential equation into first order.