pliny(t,h);

function hd = pliny(t,h);
% This function respresents the model of the Pliny's intermittent fountain
% in the form required by ODE23 and ODE45

% --- define system parameters:

R      = 0.05;          % radius of the base
r      = 0.01;          % radius of the muzzle
hlo    = 0.025;         % hight  of the orifize
hhi    = 0.1;           % hight  of the syphon
C      = 0.6;           % discharge coefficient
influx = 50e-6;         % a deciliter per second
g      = 9.81;          % gravity acceleration, m/s/s

% --- delete the following line if you are using a version of
% MATLAB prior to MATLAB 4.0
global syphon
% The global variable dry indicates whether the syphon is dry.
% If the level of water in the tank is below hlo,
% the syphon is empty.
% If the level of the water in the tank is above hhi,
% the syphon is full of water.
% If the level of the water in the tank is between hlo and hhi,
% the variable dry remain unchanged.

if h>=hhi,
  syphon=1;
end;
if h<=hlo,
  syphon=0;
end;

if syphon,
   v =C*sqrt(g*h);
   outflux = -pi*r^2*v;
else
   outflux = 0;
end;
hd =(influx+outflux)/(pi*R^2);