Error in ode arguments (line 92), error in ode45 (line 107)

In addition to the above error messages, I also get the error message 'Error in @(t,h)odefcn(t,h,A1,A2,a1,a2,K,g,V)' when I try to run the code below. Perhaps something small is amiss? Many thanks in advance.

clear all

clc

% System parameters:

A1 = 1; % (m^2) Cross-sectional area of the upper tank

A2 = 1.2; % (m^2) Cross-sectional area of the lower tank

a1 = 0.12; % (m^2) Area of the upper tank's orifice

a2 = 0.1; % (m^2) Area of the lower tank's orifice

K = 0.3; % (m^3/V*s) Pump constant

g = 9.8; % (m/s^2) Gravitational acceleration

% Initial conditions

h1_0 = 1; % (m) Initial water level in the upper tank

h2_0 = 0; % (m) Initial water level in the lower tank

V = @(t) 1 + 0.1*sin(t); % Input function (anonymous function)

%---------------------------------------------------------------

tspan = [0 50];

y0 = [h1_0 h2_0];

f = @(t,h) odefcn(t,h,A1,A2,a1,a2,K,g,V);

[t,x] = ode45(f,tspan,y0);

plot(t,x(:,1))

hold on

plot(t,x(:,2))

function dhdt = odefcn(t,h,A1,A2,a1,a2,K,g,V)

dhdt = [-(a1/A1)*sqrt(s*g*h(1)) + (K/A1)*V(t);...

(a1/A2)*sqrt(2*g*h(1)) - (a2/A2)*sqrt(2*g*h(2))];

end

Answer 1

[Sam Chak]

Hi @Thomas

There is a typo in the function odefcn().

% System parameters:
A1    = 1;      % (m^2) Cross-sectional area of the upper tank 
A2    = 1.2;    % (m^2) Cross-sectional area of the lower tank
a1    = 0.12;   % (m^2) Area of the upper tank's orifice
a2    = 0.1;    % (m^2) Area of the lower tank's orifice
K     = 0.3;    % (m^3/V*s) Pump constant
g     = 9.8;    % (m/s^2) Gravitational acceleration

% Initial conditions
h1_0  = 1;      % (m) Initial water level in the upper tank
h2_0  = 0;      % (m) Initial water level in the lower tank
V     = @(t) 1 + 0.1*sin(t); % Input function (anonymous function)

% Solving the ODE
tspan = [0 50];
y0    = [h1_0 h2_0];
f     = @(t,h) odefcn(t, h, A1, A2, a1, a2, K, g, V);
[t,x] = ode45(f, tspan, y0);

% Plotting the results
plot(t, x), grid on
legend('x_1', 'x_2')
xlabel('t'), ylabel('\bf{x}(t)')

function dhdt = odefcn(t, h, A1, A2, a1, a2, K, g, V)
    dhdt = [-(a1/A1)*sqrt(2*g*h(1)) + (K/A1)*V(t);...   % <-- s changed to 2
             (a1/A2)*sqrt(2*g*h(1)) - (a2/A2)*sqrt(2*g*h(2))];
end

[Thomas]

Hi thanks, simple typo!