Using Event Function in ode45 "Unable to perform assignment because the left and right sides have a different number of elements"

I'm trying to use an event function in this ODE. What I want this code to do, is when h=3, it switches from the first function to the second one. I'm new to matlab and this is my first time trying to use ODE45, so I'm having trouble with some of the specifics with it.

I think I have the event function set up correctly to do this, but I'm getting this error, as well as some other problems): "Unable to perform assignment because the left and right sides have a different number of elements."

(I'll provide the full error I'm getting at the bottom)

Here is my code I'm using:

global VrMax hRMax C1 Ac g C2 Hspill Lw
%% Reservoir Mass Conservation Equation (below spillway)
% dV/dt = VR-Max(h/hR-Max)^0.3 - C1Ac*sqrt(2gh)

VrMax = 1000000; % m^3
hRMax = 6; % m
C1 = 45.58; % 
Ac = .28; % m^2
g = 9.81; %m^2/s

%% Reservoir Mass Conservation Equation (over spillway)
% dV/dt = VR-Max(h/hR-Max)^0.3 - C1Ac*sqrt(2gh)+C2Lw(h-Hspill)^3/2

C2 = 0.62; 
Hspill = 3; % m
Lw = 3; %m

%% Differential Equation Parameters
t0 = 0; % Initial time of simulation 
t1 = 10; % Final time of simulation 
tspan = [t0 t1]; 

h0 = 1; % Initial condition 
h1 = 6; % Max reservoir height
hspan = [h0 h1];

options = odeset('Events',@spillway);
[t,h,t1,h1] = ode45(@masscon,tspan,hspan,options);

%% ODE Function
function dvdt = masscon(t,h)
    global VrMax hRMax C1 Ac g C2 Hspill Lw
    dvdt = zeros(2,1);
    dvdt(1) = (VrMax*(h/hRMax).^0.3) - (C1*Ac*sqrt(2*g*h))
    dvdt(2) = (VrMax*(h/hRMax).^0.3) - (C1*Ac*sqrt(2*g*h)) + (C2*Lw*(h-Hspill).^(3/2))
end
 
%% Event Function
function [position,isterminal,direction] = spillway(t,h)
position = h(3); 
isterminal = 1;  % Halt integration 
direction = 0;   % The zero can be approached from either direction
end

Thank you for any help insight to where the code is going wrong!

Full error:

Unable to perform assignment because the left and right sides have a different number of elements.

Error in final1>masscon (line 61)
    dvdt(1) = (VrMax*(h/hRMax).^0.3) - (C1*Ac*sqrt(2*g*h))

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in final1 (line 55)
[t,h,t1,h1] = ode45(@masscon,tspan,hspan,options);

[madhan ravi]

please post the latex form of the equation

Answer 1

[Walter Roberson]

h0 = 1; % Initial condition 
h1 = 6; % Max reservoir height
hspan = [h0 h1];

so hspan is a vector of length 2

[t,h,t1,h1] = ode45(@masscon,tspan,hspan,options);

hspan is being passed as initial conditions, so this ode has two state variables as input, with the conditions at any one time step being passed as the second variable of the ode function.

function dvdt = masscon(t,h)

So the second variable is being received as h. h will be a scalar of length 2 because the initial conditions hspan is length 2.

    dvdt(1) = (VrMax*(h/hRMax).^0.3) - (C1*Ac*sqrt(2*g*h))

Right hand side uses h unindexed so unless the expression manages to merge vector values somehow, the right hand side is going to be a vector. Vector values could hypothetically be merged if h was involved in an algebraic matrix multiplication or in a mrdivide or mldivide with parameters that were not scalar . For example since h will be 2 x 1, if VrMax were 1 x 2, then the * operator VrMax*expression_in_h could give a 1 x 1 result. But VrMax is a scalar, and hRMax is a scalar and so are C1 Ac and g so there is no combining done due to matrix operations. So the right hand side is calculating a vector and you are trying to store it in a scalar.