Problem with ode45 solver and initial conditions for it

I'm trying to numerically solve the optimal control problem. When I try to generalize my programm to high order systems I meet the strange problem with ode45.

When the initial conditions is written manually the ode45 correctly solves the problem. But if i taking the initial conditions from WorkSpace or as result of previous programm the solver goes wrong.

Here is the code for checking:

function PSystemSolve(psi0,x0, SysOrder,k)

%Correct working example
X0=[-1.0000  0 1.0000 0.7951];%manually entered initial condition vector 
options1=odeset('Events', @odeEvent, 'RelTol',1e-9,'AbsTol',1e-9);
[~,~,TE,XE,IE]=ode45(@NonLinSystem, [0 inf],X0, options1)

%Wrong working example
options1=odeset('Events', @odeEvent, 'RelTol',1e-9,'AbsTol',1e-9);
[~,~,TE,XE,IE]=ode45(@NonLinSystem, [0 inf],[x0 psi0], options1)

function [value, isterminal, direction]=odeEvent(t,x)
global SysOrder k
value=[x(k); x(2*SysOrder)];
isterminal=[1; 0];
direction=0;

function RPF=NonLinSystem(~,x) %The simple system
RPF=[x(2);...
    -x(2)+abs(x(4))/x(4);...
    0;...
    x(4)-x(3)];

I thought that problem is with classes of psi0 and x0, but they are double and all must be clear with solver.

If this information is not enough to answer I'll send the other part of code. I don't do this in first message, beacuse the whole program is very large.

[Walter Roberson]

Compute

[x0 psi0] - X0

and I suspect you will see some non-zero values due to round-off error in computation. For example, the 1.0000 might really be 1.0000 - 1e-16 or so.

[Ivan Khomich]

Thanks for your answer, there you are right.

But now I can't interpret this results. How the difference near 10e-6 can affect so much that problem becomes unsolvable.

The most funny thing that this simple problem at first was solved by my hand and I don't meet any computational problems, even with some round-up.

[Walter Roberson]

I think we would need the code involved in the ode, and also the exact psi0 and x0 that are needed (either attach a .m or use num2hex() to show their internal representation.)

global SysOrder k

My first assumption is always that problems with code are caused by use of global variables; I find it is not worth debugging anything else until the global variables have been eliminated.

[Ivan Khomich]

So, if it will not bother you, can you explain why global variables may cause some problems and if they tipically cause some problems in code why Matlab provide this opportunity?

Answer 1

[Ivan Khomich]

Thanks a lot, Walter.

The final problem was with the method I'd used to solve this problem, but your answer help me a lot with understanding some subtle points of how the code works in practice.