How to change a variable equation while ODE is operating

I am running a mass matrix solver using ode15s for a four-cycle engine. I have two variables, and , that are stated below:

PsiOut = sqrt(gamma/(gamma-1)*((Po/p(theta))^(2/gamma) -...
             (Po/p(theta))^((gamma+1)/gamma)));
PsiIn = sqrt(gamma/(gamma-1)*((p(theta)/Po)^(2/gamma) -...
             (p(theta)/Po)^((gamma+1)/gamma)));

p(theta) is a state variable, Po and gamma are constants. Ψ is used as a coefficient to determine how much mass is entering or leaving the system.

They are used depending on which cycle the engine is at.

At a certain point during the integration p(theta) will cross Po, which raises an issue since that causes Ψ to become imaginary.

I have two questions regarding this predicament.

1. Is the best option to solve this be to use an event function in the ode? If not, then what would be
2. If an event function is the best, how would I go about that?

I've looked at the documentation and have been working with event functions for ode, here is the code I have for:

         function [val, isterminal, dir] = PsiPositiveEventFunction(theta,p)
%%This is an event function for ode15s to interrupt the integration
%when Po = p(theta) so that Psi does not become imaginary
Po = 101324;  %Atmospheric Pressure (N/m^2)

val = p(theta) - Po;    %The value that we want to be zero for the event
isterminal = 1;     %This halts the integration if event is met
dir = 1;            %From the direction where p(theta) is growing

However, the val variable errors and it gives me back:

Array indices must be positive integers or logical values.

The bigger question, though, is how do I get back to the integration after the stoppage and continue on with the same values but with the Ψ being switched?

[Walter Roberson]

Yes, the best solution is to use an event function.

val = p(theta) - Po;    %The value that we want to be zero for the event

That line requests that you index p at location theta. That will not work.

p(theta) is a state variable

The first argument of the event function, theta, will be the current value of the independent variable, and the second argument will be the current value of the state variables. The independent variable will almost certainly have fractional components and not be suitable as an index.