Problem solving a set of 5 non-linear ODE's

i have a problem when i try to solve a set of 5 odes. i tried using ode45 and the other ode's functions and it always gives me the following error:

Warning: Failure at t=8.822299e-22. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.009266e-36) at time t. > In ode45 (line 308)

here is my code:

c=2.988*10^10; 
lambda=asin(0.9); 
M_sun=1.988*10^33;
M=M_sun*10^9; 
B=10^4; 
G=6.67*10^-8; 
r_H=G*M*(c^-2)*(1+cos(lambda)); 
Omega_H=(c^3/(2*G*M))*tan(lambda/2); 
rho_rot=(Omega_H*B)/(4*pi*c); 
rho_B=r_H; 
Gamma_thr=10^6; 
alpha=0.007297; 
e=4.8032*10^-10; 
C_1=(4/9)*alpha*Gamma_thr; 
C_2=(e*Gamma_thr^4)/(4*pi*r_H^3*rho_rot);
m_e=9.10938*10^-28; 
h_tilde_a=(Gamma_thr*m_e*c^2)/(4*pi*e*r_H^2*rho_rot); 

      function beta_a=beta_a(u_tilde_a)
      beta_a=u_tilde_a.*(Gamma_thr^-2+u_tilde_a.^2).^(-1/2);
      end     

      function Gamma_tilde_a=Gamma_tilde_a(u_tilde_a) 
      Gamma_tilde_a=(Gamma_thr^-2+u_tilde_a.^2).^(1/2); 
      end

      function E_gamma_a=E_gamma_a(u_tilde_a) 
      E_gamma_a=(2/3)*(C_2/C_1)*Gamma_tilde_a(u_tilde_a).^3;    
      end

      function alpha_tilde_a=alpha_tilde_a(u_tilde_a) 
      alpha_tilde_a=Gamma_tilde_a(u_tilde_a);
      alpha_tilde_a(alpha_tilde_a<=1)=0;
      alpha_tilde_a=alpha_tilde_a*C_1;
      end

      function q_tilde_a=q_tilde_a(n_tilde_lab_a,n_tilde_lab_plus,n_tilde_lab_minus,u_tilde_plus,u_tilde_minus) 
      q_tilde_a=(1./(2*n_tilde_lab_a)).*(n_tilde_lab_plus.*alpha_tilde_a(u_tilde_plus)+n_tilde_lab_minus.*alpha_tilde_a(u_tilde_minus));
      end

      function s_tilde_1_a=s_tilde_1_a(u_tilde_a) 
      s_tilde_1_a=C_2*(Gamma_tilde_a(u_tilde_a).^3).*u_tilde_a;  
      end

      function q_tilde_1_a=q_tilde_1_a(n_tilde_lab_a,n_tilde_lab_plus,n_tilde_lab_minus,u_tilde_plus,u_tilde_minus) 
      q_tilde_1_a=(1./(2*n_tilde_lab_a)).*(n_tilde_lab_plus.*alpha_tilde_a(u_tilde_plus).*E_gamma_a(u_tilde_plus)+n_tilde_lab_minus.*alpha_tilde_a(u_tilde_minus).*E_gamma_a(u_tilde_minus));
      end

      function D_eqs=D_eqs(x,y) 
      D_eqs=zeros(5,1);
      D_eqs(1)=y(2)-y(3)-1;
      D_eqs(2)=y(2).*Gamma_tilde_a(y(4)).*q_tilde_a(y(2),y(2),y(3),y(4),y(5)).*(1./y(4))+( 1/(Gamma_thr^2*h_tilde_a) )*( y(2)./(y(4).^2.*Gamma_tilde_a(y(4))) ).*( y(1)-s_tilde_1_a(y(4))+q_tilde_1_a(y(2),y(2),y(3),y(4),y(5))-h_tilde_a*q_tilde_a(y(2),y(2),y(3),y(4),y(5)).*y(4) );
      D_eqs(3)=y(3).*Gamma_tilde_a(y(5)).*q_tilde_a(y(3),y(2),y(3),y(4),y(5)).*(1./y(5))+( 1/(Gamma_thr^2*h_tilde_a) )*( y(3)./(y(5).^2.*Gamma_tilde_a(y(5))) ).*( -y(1)-s_tilde_1_a(y(5))+q_tilde_1_a(y(3),y(2),y(3),y(4),y(5))-h_tilde_a*q_tilde_a(y(3),y(2),y(3),y(4),y(5)).*y(5) );
      D_eqs(4)=(1/h_tilde_a)*(Gamma_tilde_a(y(4))./y(4)).*( y(1)-s_tilde_1_a(y(4))+q_tilde_1_a(y(2),y(2),y(3),y(4),y(5))-h_tilde_a*q_tilde_a(y(2),y(2),y(3),y(4),y(5)).*y(4) );
      D_eqs(5)=(1/h_tilde_a)*(Gamma_tilde_a(y(5))./y(5)).*( -y(1)-s_tilde_1_a(y(5))+q_tilde_1_a(y(3),y(2),y(3),y(4),y(5))-h_tilde_a*q_tilde_a(y(3),y(2),y(3),y(4),y(5)).*y(5) );
      end

s_tilde_initial=0; 
s_tilde_final=100;
s_tilde_step=0.005;

beta_tilde_plus_0=-10^-3; 
beta_tilde_minus_0=10^-3; 
E_tilde_par_0=1; 
n_tilde_lab_plus_0=abs(1/beta_tilde_plus_0); 
n_tilde_lab_minus_0=abs(1/beta_tilde_minus_0); 
u_tilde_plus_0=(1/Gamma_thr)*beta_tilde_plus_0*(1-beta_tilde_plus_0^2)^(-1/2); 
u_tilde_minus_0=(1/Gamma_thr)*beta_tilde_minus_0*(1-beta_tilde_minus_0^2)^(-1/2); 

[s_tilde,y_solved]=ode45(@D_eqs,[s_tilde_initial:s_tilde_step:s_tilde_final],[E_tilde_par_0 , n_tilde_lab_plus_0 , n_tilde_lab_minus_0 , u_tilde_plus_0 , u_tilde_minus_0]);

im at a loss here and clueless as to what is wrong here. i would very appreciate your help in the matter.

[Torsten]

Check the values of D_eqs(i) (i=1,...,5) returned to ODE45.

Best wishes

Torsten.

[tensorisation]

how?

[Torsten]

Insert

disp(D_eqs)

and maybe

pause(5)

at the end of function D_eqs.

Best wishes

Torsten.

[tensorisation]

ok, heres what i got from it:

1.0e+23 *

     -0.0000
      5.6675
     -5.6675
     -0.0000
     -0.0000

     1.0e+23 *

      0.0000
      5.4486
     -5.9049
     -0.0000
     -0.0000

     1.0e+23 *

      0.0000
      5.3541
     -6.0433
     -0.0000
     -0.0000

     1.0e+23 *

      0.0000
      4.9273
     -6.8841
     -0.0000
     -0.0000

     1.0e+23 *

      0.0000
      4.8560
     -7.0789
     -0.0000
     -0.0000

     1.0e+23 *

      0.0000
      4.7784
     -7.3358
     -0.0000
     -0.0000

     1.0e+23 *

      0.0000
      4.7867
     -7.3284
     -0.0000
     -0.0000

     1.0e+23 *

      0.0000
      4.3240
     -9.7815
     -0.0000
     -0.0000

     1.0e+24 *

      0.0000
      0.4166
     -1.3844
     -0.0000
     -0.0000

     1.0e+23 *

      0.0000
      3.4823
      3.0537
     -0.0000
      0.0000

     1.0e+23 *

      0.0000
      3.3319
      1.0914
     -0.0000
      0.0000

     1.0e+23 *

      0.0000
      3.2392
      1.3164
     -0.0000
      0.0000

     1.0e+24 *

      0.0000
      0.3327
     -5.5921
     -0.0000
     -0.0000

     1.0e+23 *

      0.0000
      4.6563
     -7.8431
     -0.0000
     -0.0000

     1.0e+23 *

      0.0000
      4.5975
     -8.1756
     -0.0000
     -0.0000

     1.0e+24 *

      0.0000
      0.4326
     -1.0642
     -0.0000
     -0.0000

     1.0e+24 *

      0.0000
      0.4281
     -1.1343
     -0.0000
     -0.0000

     1.0e+24 *

      0.0000
      0.4228
     -1.2426
     -0.0000
     -0.0000

     1.0e+24 *

      0.0000
      0.4231
     -1.2500
     -0.0000
     -0.0000

     1.0e+24 *

      0.0000
      0.4140
     -1.5450
     -0.0000
     -0.0000

     1.0e+24 *

      0.0000
      0.4099
     -1.8856
     -0.0000
     -0.0000

     1.0e+24 *

      0.0000
      0.3903
      1.9269
     -0.0000
      0.0000

     1.0e+23 *

      0.0000
      3.8701
      4.6120
     -0.0000
      0.0000

     1.0e+23 *

      0.0000
      3.8313
      8.1343
     -0.0000
      0.0000

     1.0e+24 *

      0.0000
      0.3832
     -1.0487
     -0.0000
     -0.0000

     1.0e+24 *

      0.0000
      0.4200
     -1.3364
     -0.0000
     -0.0000

     1.0e+24 *

      0.0000
      0.4185
     -1.3921
     -0.0000
     -0.0000

     1.0e+24 *

      0.0000
      0.4111
     -1.8010
     -0.0000
     -0.0000

     1.0e+24 *

      0.0000
      0.4098
     -1.9165
     -0.0000
     -0.0000

     1.0e+24 *

      0.0000
      0.4082
     -2.0929
     -0.0000
     -0.0000

     1.0e+24 *

      0.0000
      0.4082
     -2.1034
     -0.0000
     -0.0000

     1.0e+24 *

      0.0000
      0.4054
     -2.5751
     -0.0000
     -0.0000

     1.0e+24 *

      0.0000
      0.4041
     -3.0976
     -0.0000
     -0.0000

     1.0e+24 *

      0.0000
      0.3974
      4.1640
     -0.0000
      0.0000

     1.0e+23 *

      0.0000
      3.9626
      8.4865
     -0.0000
      0.0000

     1.0e+24 *

      0.0000
      0.3948
      1.8228
     -0.0000
      0.0000

     1.0e+24 *

      0.0000
      0.3948
     -1.4253
     -0.0000
     -0.0000

     1.0e+24 *

      0.0000
      0.4073
     -2.2445
     -0.0000
     -0.0000

     1.0e+24 *

      0.0000
      0.4068
     -2.3346
     -0.0000
     -0.0000

     1.0e+24 *

      0.0000
      0.4044
     -2.9859
     -0.0000
     -0.0000

     1.0e+24 *

      0.0000
      0.4040
     -3.1673
     -0.0000
     -0.0000

     1.0e+24 *

      0.0000
      0.4035
     -3.4392
     -0.0000
     -0.0000

     1.0e+24 *

      0.0000
      0.4035
     -3.4512
     -0.0000
     -0.0000

and so on...

the values near the end of the run:

    1.0e+32 *

      0.0000
      0.0000
     -5.1633
     -0.0000
     -0.0000

     1.0e+28 *

      0.0000
      0.0000
      3.6770
     -0.0000
      0.0000

     1.0e+28 *

      0.0000
      0.0000
      1.3345
     -0.0000
      0.0000

     1.0e+28 *

      0.0000
      0.0000
      1.4730
     -0.0000
      0.0000

     1.0e+29 *

      0.0000
      0.0000
      2.8554
     -0.0000
      0.0000

     1.0e+31 *

      0.0000
      0.0000
     -1.4894
     -0.0000
     -0.0000

     1.0e+31 *

      0.0000
      0.0000
     -1.9110
     -0.0000
     -0.0000

     1.0e+30 *

      0.0000
      0.0000
      9.9022
     -0.0000
      0.0000

     1.0e+30 *

      0.0000
      0.0000
      3.0711
     -0.0000
      0.0000

     1.0e+30 *

      0.0000
      0.0000
      4.1874
     -0.0000
      0.0000

     1.0e+31 *

      0.0000
      0.0000
     -1.7576
     -0.0000
     -0.0000

     1.0e+31 *

      0.0000
      0.0000
     -1.3687
     -0.0000
     -0.0000

     1.0e+31 *

      0.0000
      0.0000
     -1.5509
     -0.0000
     -0.0000

     1.0e+32 *

      0.0000
      0.0000
     -3.0632
     -0.0000
     -0.0000

     1.0e+30 *

      0.0000
      0.0000
     -2.4458
     -0.0000
     -0.0000

     1.0e+30 *

      0.0000
      0.0000
      2.0315
     -0.0000
      0.0000

     1.0e+29 *

      0.0000
      0.0000
      9.8965
     -0.0000
      0.0000

the high values suggest theres a problem right? but what is the cause of this and how do i fix it?

[tensorisation]

i tried running the code for a much more simplified version of the equation:

function D_eqs=D_eqs(x,y)
    D_eqs=zeros(5,1);
    D_eqs(1)=y(2)-y(3)-1;
    D_eqs(2)=( 1/(Gamma_thr^2*h_tilde_a) )*( y(2)./(y(4).^2.*Gamma_tilde_a(y(4))) ).*( y(1) );
    D_eqs(3)=( 1/(Gamma_thr^2*h_tilde_a) )*( y(3)./(y(5).^2.*Gamma_tilde_a(y(5))) ).*( -y(1) );
    D_eqs(4)=(1/h_tilde_a)*(Gamma_tilde_a(y(4))./y(4)).*( y(1) );
    D_eqs(5)=(1/h_tilde_a)*(Gamma_tilde_a(y(5))./y(5)).*( -y(1) );
    end

and i still get the same problem... so any idea whats wrong with the code?

[Walter Roberson]

It would be easier to debug if you got in the habit of using

format long g

[tensorisation]

sorry, as u can already tell by now my knowledge in matlab is very limited. can u explain in more detail what do u mean, and what exactly i need to do?

[Walter Roberson]

At the MATLAB command line, give the command

format long g

Then run your code again, and show us that output.

[tensorisation]

ok, i did that and here is what i got for the output of D_eqs:

it start with:

                          -1
            -529035124800902
             529035124800902
            264517959.177291
            264517959.177291

           -81.3803658083066
            -551187517065736
             508594554078088
             275594138.67681
            254297689.162147

           -121.716849439678
            -564110799110263
             499778815830132
            282055771.662147
            249889827.919288

           -327.953340008288
            -642593517032765
             459936381151746
            321297084.907774
             229968646.85297

           -367.508112124058
            -660782219415643
             453279178421401
            330391426.340368
            226640051.913213

           -414.478460341591
            -684760381602612
             446042661272500
            342380496.113684
            223021800.826775
and so on.. and near the end i get:

           -1415.21854942231
        1.99218494238807e+22
             374082450966633
       -9.96092471419965e+15
            187041786.612347

           -1415.21848805673
       -1.51978066686815e+22
             374082450966633
        7.59890333606427e+15
            187041786.612347

           -1415.21839554875
       -4.15523269459586e+21
             374082450966633
        2.07761634776917e+15
            187041786.612346

           -1415.21843983955
       -6.37183849829922e+21
             374082450966633
        3.18591924987222e+15
            187041786.612346

           -1415.21856181997
        1.35813146892194e+22
             374082450966633
       -6.79065734614991e+15
            187041786.612346

and by 'end' i mean the moment i get the error message:

Warning: Failure at t=9.451309e-13. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.231174e-27) at time t. > In ode45 (line 308)

[Torsten]

All we can say is: Check your equations for errors.

Best wishes

Torsten.

[tensorisation]

i found a fix for the simple case (its a case where i set the functions q_tilde_a, q_tilde_1_a to zero) by adding a constraint inside the D_eqs function.

but now if i want to solve the full problem, with q_tilde_a, q_tilde_1_a, i get an error in a certain point where there is a discontinuity in the derivatives of some of the variables because in this point the value of the functions q_tilde_a, q_tilde_1_a goes from being 0 (as it was up to this point) to some non-zero value.

so now i wanna know 2 things:

1.) if the ode functions of matlab simply cannot handle a case of discontinuity in the derivatives? is there a workaround for this?

2.) if i will simply solve up to that point , and then solve again from this point on, i will avoid the direct discontinuity in the equations and this should be a good workaround?

[Torsten]

ODE functions in general have difficulties handling discontinuities in the derivatives (not only matlab ode solvers).

2.) is the usual way to handle discontinuities.

If the time where the discontinuity occurs is known in advance, just do it as you describe.

If the time is unknown, use the "event" facility of the ODE solvers to continue your computation.

Best wishes

Torsten.

[tensorisation]

thank you for the quick reply.

"use the "event" facility of the ODE solvers" - can you explain this please?

for my case, the point isnt known in advance but what is known is that it occurs when f(x)=1 where f(x) isnt known in advance, but it is some known function of one of the functions that we solve for (let us call that solution function y(x), so what we know in advance is f(y)).

so in this case i should put something like "if(f(y)<1)" inside D_eqs and then solve again using the final values i got from the solution at f(y)<1 zone as my initial values for the zone f(y)>1?

[Torsten]

Using the event facility, you define a condition when an event occurs (in your case f(y)=1) and the measures the solver should take in this case (in your case return control to the program from which you called ODE45).

For an example, see

http://de.mathworks.com/help/matlab/math/ordinary-differential-equations.html#f1-669698

Best wishes

Torsten.

[tensorisation]

i cant really understand from thier example how to explicitly write this... i need to make sure i have the syntax correct, do i need to write:

      function event=event(y)
      event=[f(y)-1,1,0];
      end

      options=odeset('Events',@event);

and put options into the solver?

[Torsten]

Exactly.

Best wishes

Torsten.

[tensorisation]

but if i will use this method, and put for the initial conditions of the next run for the f(y)>1 zone the values i got for y at the event point, won't it just start just below the event point and then pass through the event point? or are the values automaticly given for y(x_0+eps), where x_0 is the event point, just above the point? its unclear to me.

i have another question: say i am starting to solve from 0 to x_0 then from x_0 to x_f>x_0, and i also want the solution from 0 to -x_f. i know i will have another discontinuity point roughtly at -x_0 (not exact), so do i need to solve again for 0 to -x_f as i did from 0 to x_f (while also using the event function to find the event point in the negative x side) and in the end mend the two solutions for x>0 and x<0 together? i will also need to remove the solution in x=0 from one of the solutions in order not to get the point x=0 for the mended solution twice, right?

it just that ive seen theres a fucntion called odextend that can extend my solution so instead of solving twice for both x>0 and x<0 zones, can i use instead in:

solext = odextend(sol,@D_eqs,-x_f, y_0, options)

then finally use it again to solve from the negative event point to -x_f?

[Torsten]

1. You write that the discontinuity in the dervatives comes because q_tilde_a, q_tilde_1_a goes from being 0 (as it was up to this point) to some non-zero value. So just provide q_tilde_a=q_tilde_1_a=0 before the event and q_tilde_a=q_tilde_1_a=some constant value after the event.

2. It's only senseful to use odeextend if you want to extend the solution of your ODE in the same x-direction. This is not the case in your problem. But you can easily use the ode-integrators to integrate backwards in x-direction by setting tspan appropriately. E.g. If you want to integrate from 0 to -1, just set tspan=[0 -1]. The solver will start at x=0 and end in x=-1.

Best wishes

Torsten.

[tensorisation]

1. q_tilde_a, q_tilde_1_a are not just non-zero constants after the event points, they're functions of y (a known function of course). they are defined early in my code, and they appear in the equations, so why would i have to provide them with any value? right after the event point they should simply be equal to whatever value they get for y_event (y value at the event point) , and from there on they are continuous (so theres just 1 event point in x>0 and 1 event point in x<0). anyways, it doesnt really answer what i asked.

i wanted to make sure i know how to properly get a solution from both sides of the event point. is the following method correct?: say the event point is x_event. so i solve from x=0 to some x>x_event and using the event function the integration stops at x=x_event. after this i solve again from x_event to some x=x_final, this time without using any event function. here i use the y0=y_event for the initial conditions. now related to what i asked then: is it this not gonna be problematic to integrate in the 2nd time from x=x_event to x_final? because if in the 1st integration it stopped right before the event point then now it will just continue from there and pass through the event pause, giving rise to all sorts of problems. if this is the case then how is it done properly then? and if it is not, then how does it actually skip the event point in this method?

2. ok that make sense.

[Torsten]

Usually if you use the event option, you have two phases for integration. In the first phase, you calculate with a first set of functions, parameters, etc., in the second phase with a second set. So in phase 1, you call ODE45 with a right-hand side function "D_eqs_before_event", in the second phase with a right-hand side function "D_eqs_after_event". Doing so, you can usually avoid discontinuities in the interior of the intervals of integration.

In your case, you should use a function "alpha_tilde_a_before_event" where you set

alpha_tilde_a_before_event = C1*Gamma_tilde_a(u_tilde_a)

and you should use a function "alpha_tilde_a_after_event" where you set

alpha_tilde_a_after_event = 0

although the evaluation of your function f(y) in phase 2 might indicate that you are still in phase 1 before the event.

Best wishes

Torsten.

[tensorisation]

alpha is actually non-zero after the event, and zero before the event. but i also set a condition to check for the event inside D_eqs cuz untill the event has happened in each zone x>0 and x<0 2 of the 5 y functions are meaningless and i set them to zero there:

      function D_eqs=D_eqs(x,y) 
      D_eqs=zeros(5,1);
      D_eqs(1)=y(2)-y(3)-1;  
          if x>0 
              if Gamma_tilde_a(y(5))<=1
              y(2)=0;
              D_eqs(2)=0;
              D_eqs(3)=y(3).*Gamma_tilde_a(y(5)).*q_tilde_a(y(3),y(2),y(3),y(4),y(5)).*(1./y(5))+( 1/(Gamma_thr^2*h_tilde_a) )*( y(3)./(y(5).^2.*Gamma_tilde_a(y(5))) ).*( -y(1)-s_tilde_1_a(y(5))+q_tilde_1_a(y(3),y(2),y(3),y(4),y(5))-h_tilde_a*q_tilde_a(y(3),y(2),y(3),y(4),y(5)).*y(5) );    
              y(4)=0;
              D_eqs(4)=0;
              D_eqs(5)=(1/h_tilde_a)*(Gamma_tilde_a(y(5))./y(5)).*( -y(1)-s_tilde_1_a(y(5))+q_tilde_1_a(y(3),y(2),y(3),y(4),y(5))-h_tilde_a*q_tilde_a(y(3),y(2),y(3),y(4),y(5)).*y(5) );
              end
              if Gamma_tilde_a(y(5))>1
              D_eqs(2)=y(2).*Gamma_tilde_a(y(4)).*q_tilde_a(y(2),y(2),y(3),y(4),y(5)).*(1./y(4))+( 1/(Gamma_thr^2*h_tilde_a) )*( y(2)./(y(4).^2.*Gamma_tilde_a(y(4))) ).*( y(1)-s_tilde_1_a(y(4))+q_tilde_1_a(y(2),y(2),y(3),y(4),y(5))-h_tilde_a*q_tilde_a(y(2),y(2),y(3),y(4),y(5)).*y(4) );
              D_eqs(3)=y(3).*Gamma_tilde_a(y(5)).*q_tilde_a(y(3),y(2),y(3),y(4),y(5)).*(1./y(5))+( 1/(Gamma_thr^2*h_tilde_a) )*( y(3)./(y(5).^2.*Gamma_tilde_a(y(5))) ).*( -y(1)-s_tilde_1_a(y(5))+q_tilde_1_a(y(3),y(2),y(3),y(4),y(5))-h_tilde_a*q_tilde_a(y(3),y(2),y(3),y(4),y(5)).*y(5) );
              D_eqs(4)=(1/h_tilde_a)*(Gamma_tilde_a(y(4))./y(4)).*( y(1)-s_tilde_1_a(y(4))+q_tilde_1_a(y(2),y(2),y(3),y(4),y(5))-h_tilde_a*q_tilde_a(y(2),y(2),y(3),y(4),y(5)).*y(4) );
              D_eqs(5)=(1/h_tilde_a)*(Gamma_tilde_a(y(5))./y(5)).*( -y(1)-s_tilde_1_a(y(5))+q_tilde_1_a(y(3),y(2),y(3),y(4),y(5))-h_tilde_a*q_tilde_a(y(3),y(2),y(3),y(4),y(5)).*y(5) );                       
              end           
          end 
          if x<0 
              if Gamma_tilde_a(y(4))<=1
              D_eqs(2)=y(2).*Gamma_tilde_a(y(4)).*q_tilde_a(y(2),y(2),y(3),y(4),y(5)).*(1./y(4))+( 1/(Gamma_thr^2*h_tilde_a) )*( y(2)./(y(4).^2.*Gamma_tilde_a(y(4))) ).*( y(1)-s_tilde_1_a(y(4))+q_tilde_1_a(y(2),y(2),y(3),y(4),y(5))-h_tilde_a*q_tilde_a(y(2),y(2),y(3),y(4),y(5)).*y(4) );
              y(3)=0;
              D_eqs(3)=0;
              D_eqs(4)=(1/h_tilde_a)*(Gamma_tilde_a(y(4))./y(4)).*( y(1)-s_tilde_1_a(y(4))+q_tilde_1_a(y(2),y(2),y(3),y(4),y(5))-h_tilde_a*q_tilde_a(y(2),y(2),y(3),y(4),y(5)).*y(4) );
              y(5)=0;
              D_eqs(5)=0;
              end
              if Gamma_tilde_a(y(4))>1
              D_eqs(2)=y(2).*Gamma_tilde_a(y(4)).*q_tilde_a(y(2),y(2),y(3),y(4),y(5)).*(1./y(4))+( 1/(Gamma_thr^2*h_tilde_a) )*( y(2)./(y(4).^2.*Gamma_tilde_a(y(4))) ).*( y(1)-s_tilde_1_a(y(4))+q_tilde_1_a(y(2),y(2),y(3),y(4),y(5))-h_tilde_a*q_tilde_a(y(2),y(2),y(3),y(4),y(5)).*y(4) );
              D_eqs(3)=y(3).*Gamma_tilde_a(y(5)).*q_tilde_a(y(3),y(2),y(3),y(4),y(5)).*(1./y(5))+( 1/(Gamma_thr^2*h_tilde_a) )*( y(3)./(y(5).^2.*Gamma_tilde_a(y(5))) ).*( -y(1)-s_tilde_1_a(y(5))+q_tilde_1_a(y(3),y(2),y(3),y(4),y(5))-h_tilde_a*q_tilde_a(y(3),y(2),y(3),y(4),y(5)).*y(5) );
              D_eqs(4)=(1/h_tilde_a)*(Gamma_tilde_a(y(4))./y(4)).*( y(1)-s_tilde_1_a(y(4))+q_tilde_1_a(y(2),y(2),y(3),y(4),y(5))-h_tilde_a*q_tilde_a(y(2),y(2),y(3),y(4),y(5)).*y(4) );
              D_eqs(5)=(1/h_tilde_a)*(Gamma_tilde_a(y(5))./y(5)).*( -y(1)-s_tilde_1_a(y(5))+q_tilde_1_a(y(3),y(2),y(3),y(4),y(5))-h_tilde_a*q_tilde_a(y(3),y(2),y(3),y(4),y(5)).*y(5) );    
              end
          end
      end

i finally got it working but im unsure yet if the results make sense and there arent more problems to settle. i also wonder what if i have a case where the event is repeated multiple times: i start with f(y)<1, then the event first happen when f(y)=1, then f(y) goes up and eventually down and once again we have f(y)=1 and so on for a number of times? do i need to explicitly handle each event point, or is there an efficient way for dealing with this?

[tensorisation]

now that i got it working, i ran into a strange error. i tried changing one of my parameters that appear in the equations from 1 to some higher value (say 2 for example) and i get the following error:

Error using horzcat Dimensions of matrices being concatenated are not consistent.

Error in ode15s (line 821) ieout = [ieout, ie];

Error in steady_model_cr (line 149)

im not sure what this error means, and what is the cause of it, since for exactly the same code i get no such error when one of my parameters is slightly lower.

Answer 1

[John D'Errico]

I won't try to run your code, nor will I try to figure out if there is a problem in the code, since I can't know if your code accurately reflects a realistic model. That is your job.

However, the error message indicates a classic problem - that the system is a stiff one. As importantly, the wide range of coefficients suggests this is HIGHLY probable to be a stiff system. Even the little I chose to read of your code suggests that stiff system is the VERY first thing I would consider.

Finally, I have no idea which tools you tried. Did you actually use one of the stiff solvers? If not, why not? Stiffness is the very first thing I would consider here, and ODE45 probably the last tool I'd throw at it, not the first. Note that the names of the stiff solvers in MATLAB all end in s.

[tensorisation]

like i said in one of my comments, i tried using the other ode functions like ode15s and ode23s and i still get the same error.

as i've said in another comment, i tried running the code for a much more simple version of the equations, where the equations are much simpler and i also have only 1 constant (theres 2 constants there, that are exactly equal to each other):

h_tilde_a=Gamma_thr;

      function Gamma_tilde_a=Gamma_tilde_a(u_tilde_a)       
      Gamma_tilde_a=(Gamma_thr^-2+u_tilde_a.^2).^(1/2); 
      end

      function D_eqs=D_eqs(x,y)
      D_eqs=zeros(5,1);
      D_eqs(1)=y(2)-y(3)-1;
      D_eqs(2)=( 1/(Gamma_thr^2*h_tilde_a) )*( y(2)./(y(4).^2.*Gamma_tilde_a(y(4))) ).*( y(1) );
      D_eqs(3)=( 1/(Gamma_thr^2*h_tilde_a) )*( y(3)./(y(5).^2.*Gamma_tilde_a(y(5))) ).*( -y(1) );
      D_eqs(4)=(1/h_tilde_a)*(Gamma_tilde_a(y(4))./y(4)).*( y(1) );
      D_eqs(5)=(1/h_tilde_a)*(Gamma_tilde_a(y(5))./y(5)).*( -y(1) );
      end

and i still get the same problem. so what can i possibly do to find out how to fix this problem?