close all; clear all; clc;
% dos tanques conectados, a la misma altura

% dh1/dt = F/A1-b1/A1*sqrt(h1-h2);
% dh2/dt = b1/A2*sqrt(h1-h2)-b2/A2*sqrt(h2);

% F, flujo de entrada; A1 y A2 sección transversal; b1 y b2, constantes de F = b*sqrt(h)

function hprima = dostanques(t,h,par)
  b1 = par(1);
  b2 = par(2);
  A1 = par(3);
  A2 = par(4);
  F  = par(5);
  
  hprima(1) = F/A1-b1/A1*sqrt(h(1)-h(2));
  hprima(2) = b1/A2*sqrt(h(1)-h(2))-b2/A2*sqrt(h(2));
endfunction 

b1 = 2.5; % ft^2.5/min
b2 = 5/sqrt(6); % ft^2.5/min
A1 = 5; % ft2
A2 = 10; % ft2
F  = 5; % ft3/minls
param= [b1 b2 A1 A2 F]; % vector de paramteros a usar en la función dostanques
h0 = [12 7];  % condiciones iniciales
tspan = 0:140; % vector tiempo de la corrida 

h = lsode(@(h,t) dostanques(t,h,param),h0,tspan);

plot(tspan,h)
xlabel('t (min)')
ylabel('h(ft)')