Strange error when trying to solve a set of ODEs

i ran into a strange error when trying to solve a set of odes. when my code was working fine, 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 test (line 149) (test is my file name)

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. i also noted that if assign this parameter a much higher value like 100 or 1000 the code will run without giving me error, so i have no idea what is wrong since the code was working fine. i also see no reason as to why 1 or 100 would be good values and 2 or 10 wouldn't, since this parameter (lets call it alpha) will play the following role ( this isnt exact, but close to the exact solution ): one of my function will behave as follows:

 y(x)=A*exp( alpha*(x-x_0) )+B

where A,B,x_0 are constants and alpha is my parameter. so if alpha=1,100 is fine why wouldnt alpha=2,10 be fine? makes no sense to me.

i have my suspicions as to what is causing the probelm. i opened another thread regarding things related to that, and i have not received an answer yet. so i think it would be of great help if you can answer me here:

http://www.mathworks.com/matlabcentral/answers/264266-ode-solver-event-function-problems

i am solving for x>0 and x<0 separetly because i have an initial condition for x=0, and i want the solution in both x>0 and x<0. the discontinuities are actually in the points where the function Gamma_tilde_a(y(4))=1 or Gamma_tilde_a(y(5))=1 as i specified in the event function:

    function Gamma_tilde_a=Gamma_tilde_a(u_tilde_a) % a function that expresses the normalized gamma factor (normalized with Gamma_thr) in terms of u_tilde_a %
      Gamma_tilde_a=(Gamma_thr^-2+u_tilde_a.^2).^(1/2); 
      end

      function [condition_pp,indicator_terminate_pp,direction_pp]=pair_production_and_stopping_points(x,y) % the event that specify for the ode solver when pair production occurs , applied for both the positrons and electrons after the 1st pair production had already occured %
      condition_pp=[Gamma_tilde_a(y(4))-1,Gamma_tilde_a(y(5))-1,abs(y(4))-u_tilde_a_stop_tol,abs(y(5))-u_tilde_a_stop_tol]; % for the first two components of this vector : when this value is zero, the pair production event occurs , the last two component of this vector: checking for the stopping point of the positrons (which for them will occur in the s_tilde>0 zone) and the electrons (which for them will occur in the s_tilde<0 zone) %
      indicator_terminate_pp=[0,0,1,1]; % 0 - the solver does not terminate when the event occurs , 1 - the solver does terminate when the event occurs %
      direction_pp=[0,0,0,0]; % 0 - all zeros are to be located, regardless if the event function is increasing or decreasing , +1 - only for zeros where the event function is increasing , -1 - only for zeros where the event function is decreasing %
      end

say the first event point is x_0, so in my code i first solve from x=0 to x=x_0 (using the event function to terminate the integration at x=x_0) then i solve again from x_0 to some final x value x_final. after this is done, i do the exact same thing for the x<0 zone. the error message occurs when solving for x>0 and just a little after x=x_0 (so this is before even reaching the part where it solves for x<0).

now the thing is, as i said, when i solve from x=0 to x_0 i make sure i terminate at x=x_0 using the event function, because i want to enforce a new initial condition at the point x=x_0. when i then continue to solve from x=x_0 to x_final i don't really want to terminate at any point, but the condition Gamma_tilde_a(y(4)=1 or Gamma_tilde_a(y(5))=1 may still happen so i may end up with more discontinuous points. for this reason, i used the event function and specified the events Gamma_tilde_a(y(4))=1, Gamma_tilde_a(y(5))=1 , but i didnt make it terminate. now i need to know if specifying the discontinuous points as an event in itself actually does anything in regards to dealing with the discontinuities, or do i have to terminate each time an event occurs (without knowing in advance how many times the event might occur)? so maybe i should try making a while loop that will make the integration terminate at each event point, and as long as the numel(y_event) for the y_event the ode solver returns from each time it stopped at an event point isnt zero, it will continue to integrate and terminate at each event point??

by the way, just to make things clearer, the reason that those event points are discontinuous is because i have the function:

    function alpha_tilde_a=alpha_tilde_a(u_tilde_a) % a function for alpha_tilde_a that appears in the expressions for the source terms %
        if Gamma_tilde_a(u_tilde_a)<=1
        alpha_tilde_a=0;
        else
        alpha_tilde_a=alpha_parameter;
        end
    end

which appears in my equations.

Answer 1

[Star Strider]

If it works in some situations and not others, it usually means that the row or column of the vector being added contains an empty value, so the length of the vector is less than the others. It will probably be worthwhile to use the debugging tools to see where the problem occurs.

[tensorisation]

but if this is the case, then why would it depend on the parameter value i described?

my knowledge in matlab is very limited so i have never used the debugging tools. how do i use the debugging tools to try and find out what is wrong with the code?

[Star Strider]

I don’t have your code, so I can only assume that there could be a find call or other such operation that could return an empty value. An Inf or NaN value might throw an error, but they are legitimate numerical values (as far as the IEEE standard is concerned), so I doubt they would cause that problem. Logical operations would return a logical 0 or 1 that would convert to a numerical value when used in a calculation.

For debugging help see Debugging. There are also links to demos and blogs on the Support page.

[tensorisation]

i provided the code in the comment below (in response to John)

[Star Strider]

Did you do any of the debugging steps?

What was the result?

I have no real desire to run all that when I have no idea what you’re doing. You have to take some initiative in using the debugging tools (for example dbstop).

[tensorisation]

well, the line in my code that matlab pointed to is:

   [s_tilde_positive_after_1st_pp,y_solved_for_positive_s_tilde_after_1st_pp,s_tilde_positive_pps_after_1st_pp,y_solved_for_positive_s_tilde_pps_after_1st_pp,index_positive_pps_after_1st_pp]=ode15s(@D_eqs_after_1st_pp,s_tilde_sol_span_positive_after_1st_pp,y_solved_for_positive_s_tilde_1st_pp,options_eqs_for_pps_after_1st_pp);

i dont see anything wrong with it (and as i said, for some values of alpha it doesn't give an error) so what can i do with debugging that is gonna help me find the problem?

[Star Strider]

It is extremely difficult for me to follow that code. Functions have their own workspaces, so you would have to save (or print to the Command Window) the interim values of the vector ‘D_eqs_after_1st_pp’ creates as the first step in troubleshooting your code. Perhaps just saving (or printing) the length of the vector and the (x,y) values you use to calculate it would be enough.

One problem could be that the vector you’re returning is the same name as your function, since this could create recursion problems:

D_eqs_after_1st_pp=D_eqs_after_1st_pp(x,y)

although that might not explain why it works with some parameter values and not others.

[tensorisation]

so in the end of the function:

      function D_eqs_after_1st_pp=D_eqs_after_1st_pp(x,y) % a function that describe the diffrential equations for the solver function of matlab that solves the equations, with x=s , y=[E_tilde_par ; n_tilde_lab_plus ; n_tilde_lab_minus ; u_tilde_plus ; u_tilde_minus] , D_eqs=dy/dx %
    D_eqs_after_1st_pp=zeros(5,1);
    D_eqs_after_1st_pp(1)=y(2)-y(3)-1;  
    D_eqs_after_1st_pp(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_after_1st_pp(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_after_1st_pp(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_after_1st_pp(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

ur saying i should add:

disp(y);
disp(D_eqs_before_1st_pp);

?

[Star Strider]

That is one option. I would start with:

fprintf(1, '\n\tAt x = %f\t\tLength ''D_eqs_before_1st_pp'' = %.0f\n', x, length(D_eqs_before_1st_pp));

rather than writing out the entire vector. You can re-run your code and write out the entire vector if this function turns out to be the problem. You don’t yet know what that is.

[tensorisation]

but in the begining of the function D_eqs_after_1st_pp there is:

D_eqs_after_1st_pp=zeros(5,1);

so how could it be that length(D_eqs_after_1st_pp) will be anything but 5?

by the way, the 1 in the start of the fprintf command you gave here means that it will be printed on the matlab window right? i tried putting ur line of fprintf with more digits displayed for x and here is what i got:

At x = 1.89020804506215    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020804506215    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020804506215    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020804506215    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020804506215    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020804506215    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020807326418    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020848134032    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020848134032    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020817594560    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020817594560    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020808432719    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020808432719    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020805684166    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020805684166    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020806862117    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020808040068    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020819819579    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020819819579    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020808040068    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020808040068    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020808040068    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020808040068    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020808040068    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020808040068    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020819819579    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020819819579    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020831599090    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020831599090    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020843378600    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020843378600    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020961173706    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020961173706    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020843378600    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020843378600    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020843378600    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020843378600    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020843378600    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020843378600    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020961173706    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020961173706    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020878717132    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020878717132    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020853980160    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020853980160    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020864581719    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020864581719    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020875183279    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020875183279    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020981198874    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020981198874    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020875183279    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020875183279    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020875183279    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020875183279    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020875183279    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020875183279    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020981198874    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020981198874    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020906987958    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020906987958    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020884724683    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020884724683    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020878045700    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020878045700    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020876042005    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020876042005    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020876900732    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020876900732    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020877759458    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020877759458    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020876900732    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020876900732    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020876900732    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020876900732    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020876900732    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020876900732    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020877759458    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020877759458    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020877158350    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020877415967    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020877673585    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020880249764    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020880249764    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020877673585    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020877673585    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020877673585    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020877673585    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020877673585    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020877673585    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020880249764    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020880249764    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020878446439    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020878446439    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020877905441    Length 'D_eqs_before_1st_pp' = 5

  At x = 1.89020877905441    Length 'D_eqs_before_1st_pp' = 5

i now wonder if this make sense that some x values dont change from line to line?

also, i tried printing the data into a text file by putting:

    fileID = fopen('test.txt','w');
    fprintf(fileID,'%s %s\n','s_tilde','length(D_eqs_after_1st_pp)');
    fprintf(fileID,'%f %f\n',x,length(D_eqs_after_1st_pp));
    fclose(fileID);

instead of the fprintf line, and for some reason i only got 1 line of numbers printed there. what is wrong with my usage of fprintf? (sorry for the distraction from the main topic, but my knowledge in matlab is very limited so i need to know)

aside from all this, i might have an idea as to what might cause problems in the code. it could be related to possible discontinuities in the derivatives of the function im solving for and not being able to properly treat them using the event function of the ode solver of matlab. i have already opened another thread about this issue but nobody has yet replied. if u could answer me there about the questions i presented there, it could help me find out what exactly is the problem with my code. here is a link:

http://www.mathworks.com/matlabcentral/answers/264266-ode-solver-event-function-problems

[Star Strider]

I noticed the preallocation, but had previously done this little experiment:

V = zeros(5,1);
V(3) = [];
LenV = length(V)

so if one element is for whatever reason set to ‘empty’, the length of the vector will shorten by one. That could throw the error you’re seeing.

I get the impression that your function did not throw the error while you were monitoring its activity. As for the time steps, the MATLAB ODE solvers are adaptive, so similar time steps are to be expected. (I haven’t delved into the solver codes to see exactly how they work.)

I noticed that ‘D_eqs_before_1st_pp’ has two if blocks, and no ODE solver handles significant discontinuities well. It is extremely difficult to follow your code, and while discontinuities are not differentiable, some are ‘close enough’ to not cause significant problems. If yours are significant, I would break the code up into two ode15s calls, one with x<0 and one with x>0. Use the last values of the x<0 solutions as the initial conditions for the x>0 ode15s call. I didn’t mention that before because that function did not seem to be throwing the error.

I haven’t used event functions in a while, so I can’t comment other than to say the discontinuities could very well be causing problems. I can’t say that their interactions with your event functions are causing the concatenation error.

[tensorisation]

1.) from the looks of the error matlab gave, it seems the vector that is having length problems is the one inside the ode function: Error in ode15s (line 821) ieout = [ieout, ie]; but i guess there could be many reasons for this, not just a problem with the length of one of my vectors, right?

2.) "I get the impression that your function did not throw the error while you were monitoring its activity" - what do u mean? i put the fprintf in the end of the function D_eqs_after_1st_pp, the very same function that appears in error line of my code:

Error in test (line 177) [s_tilde_positive_after_1st_pp,y_solved_for_positive_s_tilde_after_1st_pp,s_tilde_positive_pps_after_1st_pp,y_solved_for_positive_s_tilde_pps_after_1st_pp,index_positive_pps_after_1st_pp]=ode15s(@D_eqs_after_1st_pp,s_tilde_sol_span_positive_after_1st_pp,y_solved_for_positive_s_tilde_1st_pp,options_eqs_for_pps_after_1st_pp);

3.)i meant the x values themselves are the same between different lines, which means the step=0 right? that is not the same thing as having a constant step. am i missing something there? wouldnt having the x value not changing mean step=0 between the 2 lines?

4.)regarding what i asked in my previous comment, can u answer me this?: "i tried printing the data into a text file by putting:

    fileID = fopen('test.txt','w');
    fprintf(fileID,'%s %s\n','s_tilde','length(D_eqs_after_1st_pp)');
    fprintf(fileID,'%f %f\n',x,length(D_eqs_after_1st_pp));
    fclose(fileID);

instead of the fprintf line, and for some reason i only got 1 line of numbers printed there. what is wrong with my usage of fprintf?"

5.)i am solving for x>0 and x<0 separetly because i have an initial condition for x=0, and i want the solution in both x>0 and x<0. the discontinuities are actually in the points where the function Gamma_tilde_a(y(4))=1 or Gamma_tilde_a(y(5))=1 as i specified in the event function:

    function Gamma_tilde_a=Gamma_tilde_a(u_tilde_a) % a function that expresses the normalized gamma factor (normalized with Gamma_thr) in terms of u_tilde_a %
      Gamma_tilde_a=(Gamma_thr^-2+u_tilde_a.^2).^(1/2); 
      end

      function [condition_pp,indicator_terminate_pp,direction_pp]=pair_production_and_stopping_points(x,y) % the event that specify for the ode solver when pair production occurs , applied for both the positrons and electrons after the 1st pair production had already occured %
      condition_pp=[Gamma_tilde_a(y(4))-1,Gamma_tilde_a(y(5))-1,abs(y(4))-u_tilde_a_stop_tol,abs(y(5))-u_tilde_a_stop_tol]; % for the first two components of this vector : when this value is zero, the pair production event occurs , the last two component of this vector: checking for the stopping point of the positrons (which for them will occur in the s_tilde>0 zone) and the electrons (which for them will occur in the s_tilde<0 zone) %
      indicator_terminate_pp=[0,0,1,1]; % 0 - the solver does not terminate when the event occurs , 1 - the solver does terminate when the event occurs %
      direction_pp=[0,0,0,0]; % 0 - all zeros are to be located, regardless if the event function is increasing or decreasing , +1 - only for zeros where the event function is increasing , -1 - only for zeros where the event function is decreasing %
      end

say the first event point is x_0, so in my code i first solve from x=0 to x=x_0 (using the event function to terminate the integration at x=x_0) then i solve again from x_0 to some final x value x_final. after this is done, i do the exact same thing for the x<0 zone. the error message occurs when solving for x>0 and just a little after x=x_0 (so this is before even reaching the part where it solves for x<0).

now the thing is, as i said, when i solve from x=0 to x_0 i make sure i terminate at x=x_0 using the event function, because i want to enforce a new initial condition at the point x=x_0. when i then continue to solve from x=x_0 to x_final i don't really want to terminate at any point, but the condition Gamma_tilde_a(y(4)=1 or Gamma_tilde_a(y(5))=1 may still happen so i may end up with more discontinuous points. for this reason, i used the event function and specified the events Gamma_tilde_a(y(4))=1, Gamma_tilde_a(y(5))=1 , but i didnt make it terminate. now i need to know if specifying the discontinuous points as an event in itself actually does anything in regards to dealing with the discontinuities, or do i have to terminate each time an event occurs (without knowing in advance how many times the event might occur)? so maybe i should try making a while loop that will make the integration terminate at each event point, and as long as the numel(y_event) for the y_event the ode solver returns from each time it stopped at an event point isnt zero, it will continue to integrate and terminate at each event point??

by the way, just to make things clearer, the reason that those event points are discontinuous is because i have the function:

    function alpha_tilde_a=alpha_tilde_a(u_tilde_a) % a function for alpha_tilde_a that appears in the expressions for the source terms %
        if Gamma_tilde_a(u_tilde_a)<=1
        alpha_tilde_a=0;
        else
        alpha_tilde_a=alpha_parameter;
        end
    end

which appears in my equations.

[Star Strider]

This is clearly beyond my ability. I will delete my Answer in a few hours, since I’m not solving your problem.

[tensorisation]

come on seriously? and if u delete it, my comments will be deleted too? so i would have to write it all up again if by some miracle there will be some1 here willing to help me? if so, can u please not delete it? also, could u atleast answer me for 1.)-4.)? it seems like you know ur stuff atleast enough to answer those..

thanks anyways for trying, which is more than what every1 else here has done :/

[Star Strider]

OK. They would. I’ll keep it up for a week or so, in case someone else lends their expertise. You can always copy your comments (use the ‘Edit’ option to get the original text with formatting), and save them as a text file in ‘notepad’ or another basic text editor. Don’t use Word.

I’ll give a go at your questions:

1. That the problem occurs inside ode15s means that it occurs somewhere in your ODE functions. The ODE solvers attempt to append the solved values of your dependent variables to the existing ones, and since the new vector is not the same length as the previous, it throws an error. That’s my best guess, at least.

2. It may be that the function you are monitoring is not the function that is causing the error.

3. I can’t answer that. I haven’t delved into the code for the ODE solvers.

4. You would have to open the file before you called the function, and then close it afterwards. Opening it and immediately closing it would have it print only one line. I would use a .mat file instead, using save or matfile, as appropriate.

5. My impression is that an ‘event’ should interrupt the function so it can the restart with different essential parameters. You might also want to monitor the variable that is triggering the event, and see what the variables are up to that value of ‘x’.

I apologise, but that’s the best I can do. I didn’t run your code because it deals with astrophysics and cosmology that I have not a prayer of understanding with the requisite depth, so I would not be able to determine where the problem is. (I’m a physician and biomedical engineer, and while I have some experience with ODEs and the ODE solvers, need to understand the ODE itself to figure out what the problem might be. I only encounter black holes on PBS shows, and occasional reading.)

[tensorisation]

could you please just not delete it at all? it will be more convenient for me that way. i don't see what do you gain from deleting this..?

now, i have a simple but important question about the event function. i need to make sure i understand how it works. suppose my event function is:

tol=10^-9;

    function [value,indicator_terminate,direction]=my_event(x,y) 
    value=[f(y(4))-1,f(y(5))-1,abs(y(4))-tol,abs(y(5))-tol];
    indicator_terminate=[1,1,1,1];
    direction=[0,0,0,0]; 
    end

where f is some function defined earlier in the code.

i set the options for the solver as follows:

options=odeset('Events',@my_event,'Refine',20,'RelTol',10^-5,'AbsTol',10^-8);

and eventually i solve the equations using:

    [T,Y,TE,YE,IE] = solver(@odefun,tspan,y0,options)

now, since i have put:

indicator_terminate=[1,1,1,1];

inside my event function, will this mean that the solver will terminate as soon as any of the 4 event specified in:

value=[f(y(4))-1,f(y(5))-1,abs(y(4))-tol,abs(y(5))-tol];

will occur?

if so, does this means that the variable TE in:

[T,Y,TE,YE,IE] = solver(@odefun,tspan,y0,options)

will always be a vector of 1 element (which will basicly be a simple 1 value variable) ? i say this because according to my understanding (which may be wrong), if the solver wasn't set to terminate and i had several event points along the integration span, then TE would be a vector of all the T values of the event points. but could it be that TE can infact be a vector or even a matrix when using my current settings, because have more than 1 possible event that can happen at any given time (i have 4), or some other reason? something like that might explain the error that i have:

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

[Star Strider]

First, what I get out of deleting Answers is that others don’t respond to them months or years later when I’ve completely forgotten about them. I also delete nearly all my Answers that aren’t Accepted or Voted, for the same reason. Others do as well.

Second, it has been years since I’ve used event functions. I can’t even find the code I wrote that uses them so that I can test some ideas I have.

Third, unless ‘my_event’ is a nested function (function within another function), it has no idea what ‘tol’ is.

Fourth, I believe the ‘indicator_terminate’ is an ‘or’ situation, so that if any of the conditions are satisfied, the solver take the appropriate action. The documentation isn’t clear on that, and I haven’t tested it.

Fifth — and this may be important — in the documentation, all the vectors in the event functions are column vectors. You have them as row vectors, and while that may seem to be a trivial distinction, it is occasionally important in the way MATLAB evaluates them. If your events function is not working as it should, transpose the vectors to column vectors and see if that makes a difference.

[tensorisation]

1.)about your 3rd point: my entire file is an empty function: function []=test() that i run by writing my file name (test) in the matlab workspace window. so 'my_event' should be fine with 'tol' right?

2.)about your 4th point: thats what i thought it does, but it may not be true. i tried building a while loop that in each iteration will solve from the previous event point to the next event point (so it will terminate at each event point). it goes like this:

   ...
  x=[];
  y=[];

while max(x_i)<x_final 
    x_span_i=[x_event_prev,x_final] 
    [x_i,y_i,x_event_i,y_event_i,i_event_i]=ode15s(@D_eqs,x_span_i,y_event_prev,options);
    x=[x;x_i]; 
    y=[y;y_i]; 
    x_event_prev=x_event_i; 
    y_event_prev=y_event_i;
end

and i noticed that in the 2nd iteration of the loop, x_event_i is a column vector of 2 values, which are very close to each other in thier value. this obviously creates a problem (and it also trigger the error : Dimensions of matrices being concatenated are not consistent) as x_event_i is assumed to be simply the x value of the last event point when building the integration span for each iteration (using x_span_i=[x_event_prev,x_final]). so how can this be and what fix should be appropriate here?

3.)about your 5th point: can you point me to where it says it in the documentation? it just that i wasn't sure about this, so i searched about this and found:

https://www.mathworks.com/matlabcentral/newsreader/view_thread/244928

and here he actually uses row vectors, which is why i used row vectors. this isnt the trigger for my current error in the while loop, but it definitely could cause all sorts of problems, so its good you noticed that.

[Star Strider]

1. I figured as much, since you did not report that it threw an error.

2. I can’t answer that. I can’t follow your code. However if that appears to be the origin of your concatenation inconsistency, that’s where I’d look to fix it. I have no idea what the appropriate fix would be. You will likely have to go through your code and experiment. Remember, your code is both complicated and beyond my areas of expertise.

3. I simply noticed that in the documentation they use column vectors, and I know from prior experience that sometimes the orientation of the vector is important. I suggested that as a possibility. Row vectors may work, but most of the ODE functions like column vectors for everything else.

Long week, late night here.

[tensorisation]

its just a few lines of a simple while loop. in the 1st iteration of the loop i get:

x_event_i=1.89020799745197

and in the 2nd iteration the solver somehow returns the value:

x_event_i=1.89020799745198
          1.89020987224573

so now because that x_event_i is a vector, then in the next iteration of the loop, in the very first line:

x_span_i=[x_event_prev,x_final] ;

i will get an error because x_event_prev isn't a number, its a vector and the integration span isn't set correctly because its not in the form of:

x_span_i=[x_initial,x_final] ;

the thing is, i have no idea why the solver returns x_event_i as a vector. if the solver is suppost to terminate as soon as any of the 4 events i set in:

    function [value,indicator_terminate,direction]=my_event(x,y) 
    value=[f(y(4))-1;f(y(5))-1;abs(y(4))-tol;abs(y(5))-tol];
    indicator_terminate=[1;1;1;1];
    direction=[0;0;0;0]; 
    end

occurs, then the solver should return x_event_i as the location of the event where the intergration has stopped.

the only way x_event_i should be a vector is if it doesn't terminate and then it will end up returning the locations of all the events in the integration span. but in my case there should always be one event per integration zone, because the solver terminates as soon as an event occurs.

i simply don't understand why the solver returns more than 1 event point in x_event_i. if i can understand that, maybe i can move on to solving to problem, but i'm kinda stuck on that now and i'm not sure what to do.

[Star Strider]

I don’t understand it either, but one fix could be:

x_span_i=[x_event_prev(end),x_final];

That would use the very last element of ‘x_event_prev’ (if that’s what you want to do), so the ‘x_span’ vector would then have the correct dimensions. (If you want to also update ‘x_final’ each time, you could do so as well.)

[tensorisation]

1.) using this method will make the dimensions correct but won't solve the problem. maybe it is caused because somehow there are 2 event points very close to each other? maybe i need to lower the tolerance for the solver? also, say my event is f(y)=1, and f(y) changes in each iteration of the solver (because y(x) changes), how is it even guaranteed that f(y) will be precisely 1, and not skip the event point? i can have f(y)=0.99999999 then f(y)=1.00000001 and then the entire usage of event function will be rendered useless, right?

is there maybe a way for dealing with discontinuities that is better than using event function?

2.) i tried doing what you suggested using the parameters values that were 'problematic' before. i put the line

x_event_i

inside the while loop:

    ...
    x=[];
    y=[];

while max(x_i)<x_final 
    x_span_i=[x_event_prev,x_final] 
    [x_i,y_i,x_event_i,y_event_i,i_event_i]=ode15s(@D_eqs,x_span_i,y_event_prev,options);
    x=[x;x_i]; 
    y=[y;y_i]; 
    x_event_prev=x_event_i; 
    y_event_prev=y_event_i;
end

so it will print it out. what i got is that the solver finds way too many event points. say the initial event is at x_0, then right after it prints an event at x_0+0.0001 and then at x_0+0.0002 and it goes on forever printing out the event points within very small distance from each other. this not only makes the time it takes for the solver to finish be very long (it didn't even finish after 2 hours), it also doesn't make any sense that for every step there will be an event point.

3.) when i tried doing the same thing as in 2.) but without the while loop that breaks up the solution using the event points (simply solving without terminating the integration at event points) and the solver takes forever (again, it did not finish after 2 hours) so i guess it demonstrates the same behavior as in 2.).

do you have any ideas as to what might cause this sort of behavior?

[Star Strider]

My impression from your description of the event behaviour is that your ODE is for some reason getting ‘stuck’, and keeps triggering events spaced very close in terms of ‘x’. Since your code works correctly with some parameters and not with others, there is always the possibility that there is some sort of computational problem involved. My computer science background is not sufficient to address those possibilities. My knowledge of the mathematics of astrophysics and cosmology are essentially nonexistent.

I am long since out of all ideas as to what could be causing your problem, and have long since exhausted my expertise. I was hoping that one of the other contributors with significant experience and expertise in the differential equation solvers would chime in with other ideas, but no one has.

[tensorisation]

well, the cause for trouble im getting when trying to solve for different parameters must be the discontinuities. the problem is that it seems that simply using an event function to terminate the integration and then restart it isn't working, so i wonder if there is a better way to handle discontinuities? for a problem not that uncommon i would've expected to find a solution by now :/.

the best i can do now is to try and replace the discontinuous function with some continuous one, in order to get some results. but this is hardly solves the problem, because that is not what i need.

[Star Strider]

The only way to handle discontinuities is to avoid them. I believe I previously suggested integrating to a discontinuity, stopping, and restarting with the previous values (or slightly changed values for them) as the new initial conditions. That is how the event functions are supposed to work.

I doubt any ODE solver handles discontinuities well, and the MATLAB solvers, for all their robustness, do not.

[tensorisation]

yeah that is precisely what i have done in my code in the while loop i showed before, but it doesn't seem to work properly (like i said, when i change my parameters i get the problems i mentioned earlier).

maybe my way of trying to skip the discontinuities using the event function is wrong?

for me, the discontinuities = my events , so say it happens for f(y)=1 i set the event f(y)=1 to terminate the integration, and then i start integrating again from the point i've stopped (using x_event as the starting point, and y_event as the initial condition). this process happens at each iteration of the while loop i showed earlier.

in addition, as i've already mentioned before, maybe it is possible that starting the next integration right at the point where the previous integration terminated (at the event point) and not slightly after, could cause problems? this is due to the fact that its unclear if the integration was terminated and returned x_event , y_event values that are SLIGHTLY AFTER the event point, as i would want, or SLIGHTLY BEFORE the event and in that case wouldn't i bump into the same event point all over again? in the latter case, it isn't really clear what would be the step i should add to x_event in order to get the most accurate results.

if you had an unknown number of discontinuities points at f(y)=1, how would deal with them? im asking this kind of out of desperation because as i've tried according to all i know using the while loop i showed before and it didn't work, but maybe im completely missing something..

[Star Strider]

I quite definitely do not understand the physics you are modeling, so I can only give a qualitative reply. I don’t know what the dynamics of your system are, or what your discontinuities represent. You could start after each discontinuity by changing specific values of the previous solution that characterise the discontinuity to different values, or starting over with completely new initial conditions (not related to the previous solution) to see where the dynamics of your system lead.

For an ‘unknown number of discontinuities at f(y)=1’, my first approach would be to determine the number of discontinuities, and the variables they represented. I would want to know what was going on in the system I was modeling. (For example, would taking the eigenvalues of the system in the region of the discontinuities be appropriate?) Then I would do what I could to develop a way of dealing with them.

If you are doing original research with your system (rather than reproducing previously published results), it could be that your model does not correctly describe the actual system, or the specific parameters you’re using that do not work (as opposed to the others that do) are not realistic or appropriate.

[tensorisation]

my problem is solely related to numerics/matlab programming, so in essence one does not need to know a thing about what my code is about in order to provide an answer for my problem: how to deal with discontinuities when integrating an ODE that has discontinuities in matlab?

when i tried using event function to stop the integration and then restart it, i couldn't get it to work properly, so either i didn't do it correctly (if so, i would like to know what is wrong with what i have done so i can fix it- i've already described in details what i did), or there is a different way to handle this (in this case i would like to know what is it i need to do).

the reason i said 'unknown number of discontinuities' is because the the solution y(x) obviously depens on the values of the parameters i use, and it is obviously gonna affect the event function f(y). so just as i do not know in advance the location of the event points (because i do not know in advance y(x), and the event depens on y), then i also do not know in advance the number of events that would occur. this information isn't really what is missing and is not the focus of the problem i am dealing with, because i am using the event function to locate the discontinuities as i integrate, so all of this shouldn't bother me.

in addition, it doesn't really matter if my model is physical or not. everything that is needed to be known in order to provide a solution to my numeric/matlab programming problem is that a i have a set of ODE with discontinuities. it doesn't matter how i derived the equations or what physical system they represent.

at this point, im frustrated that its been several weeks and i am still stuck on this problem. do you think it's a good idea to open a new thread about this?

[Star Strider]

‘Do you think it's a good idea to open a new thread about this?’

Yes. This has evolved over time to something far different from your original Question, and it’s quite likely others haven’t looked at it.

Be as specific as you can in your Question, and attach all your code.

[tensorisation]

my new thread is open for several days now and so far i havn't gotten a single reply. it looks like it's not gonna change any time soon. are there any sort of moderators that are responsible for providing answers? if so, i would like to know how to contact them so i can try and ask thier assistance. if not, what do you suggest i do?

[Star Strider]

That’s the problem with difficult Questions. There are a few MathWorks staffers who show up here, but the rest of us are unpaid volunteers.

You can use the ‘Contact Support’ link (See ‘Contact Us’ at the top right corner of this page), and see if MathWorks can help. They will if it’s a question about a specific function (in this instance the ‘event’ option in your ODE solver, if that’s actually the problem), so putting it in those terms could help. Reference the URL to this thread in your question to them so they don’t have to search for it and so you don’t have to repeat all of it. I haven’t seen your other post, so you may want to reference it as well.

That’s the best I can do.

[Nicolas Mira Gebauer]

Hi, I am getting a very similar problem, although 4 and a half years later... How did you solve your problem?

Thank you very much

Answer 2

[John BG]

could you please start by hanging the complete code you use?

the one that results into the error code you have already described in detail.

Awaiting answer

John

[e-mail]

[tensorisation]

well, in my experience if i give my complete code people don't bother reading it, but here you go:

format long g; % using a display format that will give an output in a more convenient manner, in case there is a check/test that we want to make %
% natrual constants %
c=2.988*10^10; % speed of light in cm/s %
e=4.8032*10^-10; % the electrical charge of the electron, in units of statC %
m_e=9.10938*10^-28; % the mass of the electron, in units of g %
hbar=1.054572*10^-27; % the reduced Plank constant, in units of erg*s %
alpha=e^2/(hbar*c); % the fine structure constant, a unitless quantity %
G=6.67*10^-8; % Newton's gravitational constant, in units of dyne*cm^2*g^-2 %
M_sun=1.988*10^33; % the mass of the sun in g %
% parameters/constants related to the blackhole itself %
lambda=asin(0.9); % sin(lambda) represents the normalized angular momentum of the blackhole %
M=M_sun*10^9; % the mass of the blackhole %
r_H=G*M*(c^-2)*(1+cos(lambda)); % the radius for the blackhole horizon, in units of cm %
Omega_H=(c^3/(2*G*M))*tan(lambda/2); % the angular velocity of the blackhole horizon, in units of s^-1 %
% additional parameters/constants %
B=10^4; % the typical magnetic field in our area of interest, in units of G (Gauss)%
rho_rot=(Omega_H*B)/(4*pi*c); % the Goldreich-Julian charge density, in units of statC/cm^3 %
rho_B=r_H; % the curvature radius of the magnetic field lines, in units of cm %
omega_B=sqrt((4*pi*e*rho_rot)/m_e); % a characteristic frequency scale in our model %
L_B=c/omega_B; % a characteristic length scale in our model %
Gamma_thr=10^6; % the effective Lorentz factor threshold for the pair production process %
C_1=(4/9)*alpha*L_B*Gamma_thr*r_H^-1; % one of the dimensionless constants i introduced in the expressions for the source terms %
C_2=(e*Gamma_thr^4)/(4*pi*r_H^2*L_B*rho_rot); % one of the dimensionless constants i introduced in the expressions for the source terms %
h_tilde_a=Gamma_thr; % the normalized enthalpy per particle for each species of fluid (in the rest frame of the fluid) %

      function beta_a=beta_a(u_tilde_a) % a function that expresses the normalized velocity (normalized with c) in terms of 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) % a function that expresses the normalized gamma factor (normalized with Gamma_thr) in terms of 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) % a function for the typical energy of a CR (curvature radiation, which is also called synchrotron radiation) photon in terms of u_tilde_a %
      E_gamma_a=(2/3)*(C_2/C_1)*Gamma_tilde_a(u_tilde_a).^3;    
      end

% changing these 2 parameters in order to find a steady state solution %
alpha_parameter=1;
xi=10^-2;

      function alpha_tilde_a=alpha_tilde_a(u_tilde_a) % a function for alpha_tilde_a that appears in the expressions for the source terms %
          if Gamma_tilde_a(u_tilde_a)<=1
          alpha_tilde_a=0;
          else
          alpha_tilde_a=alpha_parameter;
          end
      end
  %}

%{
    function alpha_tilde_a=alpha_tilde_a(u_tilde_a) % a function for alpha_tilde_a that appears in the expressions for the source terms %
        if Gamma_tilde_a(u_tilde_a)<=1
        alpha_tilde_a=0;
        else
        alpha_tilde_a=C_1*Gamma_tilde_a(u_tilde_a);
        end
    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) % a function for the source term q_tilde_a %
    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) % a function for the source term s_tilde_1_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) % a function for the source term q_tilde_1_a %
      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_before_1st_pp=D_eqs_before_1st_pp(x,y) % a function that describe the diffrential equations for the solver function of matlab that solves the equations, with x=s , y=[E_tilde_par ; n_tilde_lab_plus ; n_tilde_lab_minus ; u_tilde_plus ; u_tilde_minus] , D_eqs=dy/dx %
      D_eqs_before_1st_pp=zeros(5,1);
          if x>0 % for s_tilde>0 where there are no initial positrons, we use the constraint that: n_tilde_plus=0 , u_tilde_plus=0 as long as Gamma_tilde_minus<=1 (as long as there is no pair production) %
              y(2)=0;
              y(4)=0;
              D_eqs_before_1st_pp(1)=y(2)-y(3)-1; 
              D_eqs_before_1st_pp(2)=0;
              D_eqs_before_1st_pp(3)=( 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)) );    
              D_eqs_before_1st_pp(4)=0;
              D_eqs_before_1st_pp(5)=(1/h_tilde_a)*(Gamma_tilde_a(y(5))./y(5)).*( -y(1)-s_tilde_1_a(y(5)) );   
          end 
          if x<0 % for s_tilde<0 where there are no initial electrons, we use the constraint that: n_tilde_minus=0 , u_tilde_minus=0 as long as Gamma_tilde_plus<=1 (as long as there is no pair production) %
              y(3)=0;
              y(5)=0;
              D_eqs_before_1st_pp(1)=y(2)-y(3)-1;             
              D_eqs_before_1st_pp(2)=( 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)) );
              D_eqs_before_1st_pp(3)=0;
              D_eqs_before_1st_pp(4)=(1/h_tilde_a)*(Gamma_tilde_a(y(4))./y(4)).*( y(1)-s_tilde_1_a(y(4)) );  
              D_eqs_before_1st_pp(5)=0;
          end
      end

      function D_eqs_after_1st_pp=D_eqs_after_1st_pp(x,y) % a function that describe the diffrential equations for the solver function of matlab that solves the equations, with x=s , y=[E_tilde_par ; n_tilde_lab_plus ; n_tilde_lab_minus ; u_tilde_plus ; u_tilde_minus] , D_eqs=dy/dx %
      D_eqs_after_1st_pp=zeros(5,1);
      D_eqs_after_1st_pp(1)=y(2)-y(3)-1;  
      D_eqs_after_1st_pp(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_after_1st_pp(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_after_1st_pp(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_after_1st_pp(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; % the value of s_tilde from which we start the solution. on this value of s_tilde we have the boundary conditions. this is the injection point. %
s_tilde_final=50; % the value we choose for the final s_tilde %

% the boundary conditions at s_tilde=s_tilde_initial %
E_tilde_par_max=(e*L_B*B)/(m_e*c^2); % the maximum value for E_tilde_par obtained form setting E_par=B %
%xi=10^-2; % a parameter for characterizing E_tilde_par_0 relative to E_tilde_par_max %
E_tilde_par_0=-xi*E_tilde_par_max; % this is E_tilde_par(s_tilde=s_tilde_initial) %

Gamma_tilde_plus_0=10/Gamma_thr; % the value we choose for Gamma_tilde_plus at s_tilde=0 %
Gamma_tilde_minus_0=10/Gamma_thr; % the value we choose for Gamma_tilde_minus at s_tilde=0 %
beta_tilde_plus_0=-(1-Gamma_thr^-2*Gamma_tilde_plus_0^-2)^(1/2); % the value for beta_tilde_plus at s_tilde=0 %
beta_tilde_minus_0=(1-Gamma_thr^-2*Gamma_tilde_minus_0^-2)^(1/2); % the value for beta_tilde_minus at s_tilde=0 %
u_tilde_plus_0=Gamma_tilde_plus_0*beta_tilde_plus_0; % the value for u_tilde_plus at s_tilde=0 %
u_tilde_minus_0=Gamma_tilde_minus_0*beta_tilde_minus_0; % the value for u_tilde_plus at s_tilde=0 %

N_tilde_0_plus=-10^-2; % The normalized particle flux at s_tilde=0 for the positrons %
N_tilde_0_minus=10^-2; % The normalized particle flux at s_tilde=0 for the electrons %
n_tilde_lab_plus_0=N_tilde_0_plus/beta_tilde_plus_0; % this is n_tilde_lab_plus(s_tilde=s_tilde_initial) %
n_tilde_lab_minus_0=N_tilde_0_minus/beta_tilde_minus_0; % this is n_tilde_lab_minus(s_tilde=s_tilde_initial) %

y_0=[E_tilde_par_0 , n_tilde_lab_plus_0 , n_tilde_lab_minus_0 , u_tilde_plus_0 , u_tilde_minus_0]; % putting all the initial conditions together as a vector for the solver %
s_tilde_sol_span_positive=[s_tilde_initial,s_tilde_final]; % solution range for s_tilde>0 %
s_tilde_sol_span_negative=[s_tilde_initial,(-s_tilde_final)]; % solution range for s_tilde<0 % 

      function [condition_1st_pp_for_positive_s_tilde,indicator_terminate_1st_pp_for_positive_s_tilde,direction_1st_pp_for_positive_s_tilde]=pair_production_1st_for_positive_s_tilde(x,y) % the event that specify for the ode solver when the 1st pair production occurs , applied for the electrons in the s_tilde>0 zone %
      condition_1st_pp_for_positive_s_tilde=Gamma_tilde_a(y(5))-1; % when this value is zero, the pair production event occurs %
      indicator_terminate_1st_pp_for_positive_s_tilde=1; % for this value the solver terminate when the event occurs %
      direction_1st_pp_for_positive_s_tilde=1; % all zeros are to be located when the event function is increasing %
      end

      function [condition_1st_pp_for_negative_s_tilde,indicator_terminate_1st_pp_for_negative_s_tilde,direction_1st_pp_for_negative_s_tilde]=pair_production_1st_for_negative_s_tilde(x,y) % the event that specify for the ode solver when the 1st pair production occurs , applied for the positrons in the s_tilde<0 zone %
      condition_1st_pp_for_negative_s_tilde=Gamma_tilde_a(y(4))-1; % when this value is zero, the pair production event occurs %
      indicator_terminate_1st_pp_for_negative_s_tilde=1; % for this value the solver terminate when the event occurs %
      direction_1st_pp_for_negative_s_tilde=1; % all zeros are to be located when the event function is increasing %
      end

u_tilde_a_stop_tol=10^-9; % setting a tolerance for u_tilde_a from which below the particles will be considered to have reached a complete stop, triggering the termination of the code run (this is for preventing possible problems that might arise from not terminating the code run) %

      function [condition_pp,indicator_terminate_pp,direction_pp]=pair_production_and_stopping_points(x,y) % the event that specify for the ode solver when pair production occurs , applied for both the positrons and electrons after the 1st pair production had already occured %
      condition_pp=[Gamma_tilde_a(y(4))-1,Gamma_tilde_a(y(5))-1,abs(y(4))-u_tilde_a_stop_tol,abs(y(5))-u_tilde_a_stop_tol]; % for the first two components of this vector : when this value is zero, the pair production event occurs , the last two component of this vector: checking for the stopping point of the positrons (which for them will occur in the s_tilde>0 zone) and the electrons (which for them will occur in the s_tilde<0 zone) %
      indicator_terminate_pp=[0,0,1,1]; % 0 - the solver does not terminate when the event occurs , 1 - the solver does terminate when the event occurs %
      direction_pp=[0,0,0,0]; % 0 - all zeros are to be located, regardless if the event function is increasing or decreasing , +1 - only for zeros where the event function is increasing , -1 - only for zeros where the event function is decreasing %
      end

options_eqs_for_positive_s_tilde_1st_pp=odeset('Events',@pair_production_1st_for_positive_s_tilde,'Refine',20,'RelTol',10^-5,'AbsTol',10^-8); % special options for the ode solver , applied for the s_tilde>0 zone before the electrons cause the 1st pair production event %
options_eqs_for_negative_s_tilde_1st_pp=odeset('Events',@pair_production_1st_for_negative_s_tilde,'Refine',20,'RelTol',10^-5,'AbsTol',10^-8); % special options for the ode solver , applied for the s_tilde<0 zone before the positrons cause the 1st pair production event %
options_eqs_for_pps_after_1st_pp=odeset('Events',@pair_production_and_stopping_points,'Refine',20,'RelTol',10^-5,'AbsTol',10^-8); % special options for the ode solver , applied for both the positrons and electrons after the 1st pair production had already occured %

% solving the equations to get E_tilde_par(s_tilde) , n_tilde_lab_plus(s_tilde) , n_tilde_lab_minus(s_tilde) , u_tilde_plus(s_tilde) , u_tilde_minus(s_tilde) %

% obtaining the solution for s_tilde>0 %    
[s_tilde_positive_before_1st_pp,y_solved_for_positive_s_tilde_before_1st_pp,s_tilde_positive_1st_pp,y_solved_for_positive_s_tilde_1st_pp,index_positive_1st_pp]=ode15s(@D_eqs_before_1st_pp,s_tilde_sol_span_positive,y_0,options_eqs_for_positive_s_tilde_1st_pp); % solving from s_tilde_initial=0 up to when Gamma_tilde_minus=1 (where the 1st pair production process occurs) %
numel(y_solved_for_positive_s_tilde_1st_pp) % trying to see if theres a problem here %
if numel(y_solved_for_positive_s_tilde_1st_pp)>0 % checking to see that the solver stopped at the event point in the s_tilde>0 zone %
    s_tilde_sol_span_positive_after_1st_pp=[s_tilde_positive_1st_pp,s_tilde_final]; % solution range, from the point of the 1st pair prodution process up to s_tilde_final, at the s_tilde>0 zone %
    y_solved_for_positive_s_tilde_1st_pp=[y_solved_for_positive_s_tilde_1st_pp(1) , y_solved_for_positive_s_tilde_1st_pp(2) , y_solved_for_positive_s_tilde_1st_pp(3) , y_solved_for_positive_s_tilde_1st_pp(5) , y_solved_for_positive_s_tilde_1st_pp(5)]; % at the point of the 1st pair production process in the s_tilde>0 zone, we assign the condition u_tilde_plus=u_tilde_minus %
    [s_tilde_positive_after_1st_pp,y_solved_for_positive_s_tilde_after_1st_pp,s_tilde_positive_pps_after_1st_pp,y_solved_for_positive_s_tilde_pps_after_1st_pp,index_positive_pps_after_1st_pp]=ode15s(@D_eqs_after_1st_pp,s_tilde_sol_span_positive_after_1st_pp,y_solved_for_positive_s_tilde_1st_pp,options_eqs_for_pps_after_1st_pp); % solving from the point of the 1st pair prodution process up to s_tilde_final, at the s_tilde>0 zone %
end

% obtaining the solution for s_tilde<0 %
[s_tilde_negative_before_1st_pp,y_solved_for_negative_s_tilde_before_1st_pp,s_tilde_negative_1st_pp,y_solved_for_negative_s_tilde_1st_pp,index_negative_1st_pp]=ode15s(@D_eqs_before_1st_pp,s_tilde_sol_span_negative,y_0,options_eqs_for_negative_s_tilde_1st_pp); % solving from s_tilde_initial=0 up to when Gamma_tilde_plus=1 (where the 1st pair production process occurs) %
numel(y_solved_for_negative_s_tilde_1st_pp) % trying to see if theres a problem here %
if numel(y_solved_for_negative_s_tilde_1st_pp)>0 % checking to see that the solver stopped at the event point in the s_tilde<0 zone %
    s_tilde_sol_span_negative_after_1st_pp=[s_tilde_negative_1st_pp,(-s_tilde_final)]; % solution range, from the point of pair prodution down to -s_tilde_final, at the s_tilde<0 zone %
    y_solved_for_negative_s_tilde_1st_pp=[y_solved_for_negative_s_tilde_1st_pp(1) , y_solved_for_negative_s_tilde_1st_pp(2) , y_solved_for_negative_s_tilde_1st_pp(3) , y_solved_for_negative_s_tilde_1st_pp(4) , y_solved_for_negative_s_tilde_1st_pp(4)]; % at the point of the 1st pair production process in the s_tilde<0 zone, we assign the condition u_tilde_minus=u_tilde_plus %
    [s_tilde_negative_after_1st_pp,y_solved_for_negative_s_tilde_after_1st_pp,s_tilde_negative_pps_after_1st_pp,y_solved_for_negative_s_tilde_pps_after_1st_pp,index_negative_pps_after_1st_pp]=ode15s(@D_eqs_after_1st_pp,s_tilde_sol_span_negative_after_1st_pp,y_solved_for_negative_s_tilde_1st_pp,options_eqs_for_pps_after_1st_pp); % solving from the point of the 1st pair prodution process down to -s_tilde_final, at the s_tilde<0 zone %
end

if numel(y_solved_for_positive_s_tilde_1st_pp)>0 % checking to see that the solver stopped at the event point in the s_tilde>0 zone %
s_tilde_positive=[s_tilde_positive_before_1st_pp ; s_tilde_positive_after_1st_pp]; % mending the solution for before and after the 1st pair production event, at the s_tilde>0 zone %
y_solved_for_positive_s_tilde=[y_solved_for_positive_s_tilde_before_1st_pp ; y_solved_for_positive_s_tilde_after_1st_pp];
else
s_tilde_positive=s_tilde_positive_before_1st_pp;
y_solved_for_positive_s_tilde=y_solved_for_positive_s_tilde_before_1st_pp;
end

if numel(y_solved_for_negative_s_tilde_1st_pp)>0 % checking to see that the solver stopped at the event point in the s_tilde<0 zone %
s_tilde_negative=[s_tilde_negative_before_1st_pp ; s_tilde_negative_after_1st_pp]; % mending the solution for before and after the 1st pair production event, at the s_tilde<0 zone %
y_solved_for_negative_s_tilde=[y_solved_for_negative_s_tilde_before_1st_pp ; y_solved_for_negative_s_tilde_after_1st_pp];
else
s_tilde_negative=s_tilde_negative_before_1st_pp;
y_solved_for_negative_s_tilde=y_solved_for_negative_s_tilde_before_1st_pp;
end

% mending the solution for s_tilde>0 and s_tilde<0 : getting the values for E_tilde_par(s_tilde) , n_tilde_lab_plus(s_tilde) , n_tilde_lab_minus(s_tilde) , u_tilde_plus(s_tilde) , u_tilde_minus(s_tilde) made up from the solution in s_tilde>0 and in s_tilde<0 , while also not taking the values at s_tilde=0 twice %   
s_tilde=[ flipud(s_tilde_negative(2:length(s_tilde_negative))) ; s_tilde_positive ];
E_tilde_par=[ flipud( y_solved_for_negative_s_tilde(2:length(y_solved_for_negative_s_tilde(:,1)),1) ) ; y_solved_for_positive_s_tilde(:,1) ];
n_tilde_lab_plus=[ flipud( y_solved_for_negative_s_tilde(2:length(y_solved_for_negative_s_tilde(:,2)),2) ) ; y_solved_for_positive_s_tilde(:,2) ];
n_tilde_lab_minus=[ flipud( y_solved_for_negative_s_tilde(2:length(y_solved_for_negative_s_tilde(:,3)),3) ) ; y_solved_for_positive_s_tilde(:,3) ];
u_tilde_plus=[ flipud( y_solved_for_negative_s_tilde(2:length(y_solved_for_negative_s_tilde(:,4)),4) ) ; y_solved_for_positive_s_tilde(:,4) ];
u_tilde_minus=[ flipud( y_solved_for_negative_s_tilde(2:length(y_solved_for_negative_s_tilde(:,5)),5) ) ; y_solved_for_positive_s_tilde(:,5) ];

[John BG]

If you don't supply THE code that brought you to the error you are seeking help for, the readers cannot reproduce the problem. In the script you supplied,

1.- You are defining functions with function names exactly the same as the output variables of such functions. Do you really expect such habit to have your script, or any script, working without errors?

2.- The code you supplied give basic syntax errors. removed some of them, gave up 15min later.

3.- The following are variables NOT defined in the supplied script:

h_tilde_a=Gamma_thr
u_tilde_a=Gamma_thr
n_tilde_lab_a=Gamma_thr
n_tilde_lab_plus=Gamma_thr
u_tilde_plus= Gamma_thr
n_tilde_lab_minus= Gamma_thr

in an attempt to have the supplied script reaching the error originating this question I gave the undefined variables the above errors, no way through the syntax errors.

4.- Something as basic as the variable x, not defined I did x=xi, and at this point is a good moment to ask, noch mal:

please hang the code that produces the error originating this question, not a chunk of it, not an earlier version that produces earlier errors, not the code without variables ..

5.- Would it be possible to make the script available as attached file instead?

And for the record: MATLAB CENTRAL readers, or for the case, any reader of any document, ever, tend to give up reading when either

1.- do not want to read further: not interested, boring, already know, something better to read

2.- find it difficult to read further.

my guess your question falls in case 2.

Regards

John

[e-mail]