MATLAB: Event function with odesolver error

differential equationsMATLAB

Hello,
I am trying to integrate a vector E subject to a first order differential equation. I want the integration to stop when a certain threshold is reached. Therefore I have constructed the following functions:
d=@(s,E)dPdEfull(s,E,t,U0,H0,mu,initial,final)
s is my integration variable. E is a vector in time, with a predetermined number of components. My event function is
e=@(s,E)Pifevent(s,E,t,U0,H0,mu)
with Pifevent defining the value, isterminal, and direction. Now when I evaluate both functions separately with some E vector, they work fine. But when I try to solve by using
[s E]=ode45(d,[0 10],E0,options);
with
options=odeset('Events',e);
and E0 the initial E vector, I get an error message
Attempted to access E(2); index out of bounds because numel(E)=1.
I don't know why this happens, especially since when I try to integrate without the event function and just do
    [s E]=ode45(d,[0 10],E0);
it works fine. It seems that the event function is not reading E as a vector. Maybe it is trying to evaluate the event function before the first integration supplies a vector E? I don't know how to fix this.
Any help would be appreciated.
Edit: I'm posting more code as per Matt's request.
d(0,E0)
returns a column vector of size 101×1, I won't post the output since it's too long.
[a,b,c]=e(0,E0)
returns a=-.2923, b=1, and c=0. Does that help?
Edit 2:E0 is 1×101. The code for dPdEfull isn't that bad, I'll go ahead and post it and Pifevent:
function dEds = dPdEfull(s, E, t, U0, H0, mu, initial, final)
dEds=zeros([1 numel(E)]);
for l=1:length(t)
  m=min(t):t(2)-t(1):t(l);
  a=Schrodinger(t,U0,H0,mu,E);
  c=Schrodinger(m,U0,H0,mu,E);
  b=c';
  x=a(final,initial);
  Y=a*b*mu*c;
  y=Y(final,initial);
  dEds(l)=2*imag(x*y);
end
dEds=dEds(:);
end
and Pifevent:
function [value, isterminal, direction] = Pifevent(s,E,t,U0,H0,mu)
u=Schrodinger(t,U0,H0,mu,E(end,:));
v=u(2,3);
v=abs(v);
v=v^2;
value=v-0.99;
isterminal=1;
direction=0;
end
The function Schrodinger just takes in those arguments and returns a unitary matrix (3×3 in this example). By the way t is also a vector.

Best Answer

Ah. I think I have an idea now. Check the error message trace. I'm guessing the error actually occurs within Schrodinger. When ode45 passes E into the ODE function dPdEfull or the event function Pifevent, it passes just the current timestep's value, so E is a column vector (even though the E returned by ode45 is in the form of a matrix). Inside dPdEfull you call Schrodinger:
a=Schrodinger(t,U0,H0,mu,E);
and here E is a vector and all is well. But inside Pifevent you call
u=Schrodinger(t,U0,H0,mu,E(end,:));
Because E is a (column) vector, E(end,:) is the last row, which is a scalar. If Schrodinger expects a vector as the last argument, that could cause an error like you're getting.
Related Question