function xdot = odequadtanklin(t,x,flag,parvec,u,d)
%
% b.w. bequette - 18 Nov 09
% linear quadruple tank problem
% if parvec(1) == 1, then operating point 1
% disturbances are perturbations to manipulated inputs
% 
% parameters (continuous time, state space)
% 
  A1 = 28;    A3 = A1; % cross-sectional area
  A2 = 32;    A4 = A2; % cross-sectional area
  a1 = 0.071; a3 = a1;
  a2 = 0.057; a4 = a2;
  beta = 0.5;
  grav = 981;
  if parvec(1) == 1;
% parameters for operating point 1 (minimum phase)
    k1 = 3.33;  k2 = 3.35;
    gamma1 = 0.7;  gamma2 = 0.6;
    h1s = 12.4; h2s = 12.7;
    h3s = 1.8;  h4s = 1.4;
    v1s = 3;    v2s = 3;
  else
% operating point 2 (nonminimum phase) -- these are the first simulations
    k1 = 3.14;  k2 = 3.29;
    gamma1 = 0.43;  gamma2 = 0.34;
    h1s = 12.6; h2s = 13.0;
    h3s = 4.8;  h4s = 4.9;
    v1s = 3.15;    v2s = 3.15;
  end
% time constants
  tau1 = (A1/a1)*sqrt(2*h1s/grav);
  tau2 = (A2/a2)*sqrt(2*h2s/grav);
  tau3 = (A3/a3)*sqrt(2*h3s/grav);
  tau4 = (A4/a4)*sqrt(2*h4s/grav);
% continuous state space model
  aplant = [-1/tau1,0,A3/(A1*tau3),0;0,-1/tau2,0,A4/(A2*tau4);0,0,-1/tau3,0;0,0,0,-1/tau4];
  bplant = [gamma1*k1/A1,0;0,gamma2*k2/A2;0,(1-gamma2)*k2/A3;(1-gamma1)*k1/A4,0];
%  the states are the tank heights, in deviation variables
%
  xdot = aplant*x + bplant*u + bplant*d;
  