%
%   Tanks2_unforced.m    Numerical Integration of System of ODEs (IVP)
%
%   Sample problem from in-class notes -- homogeneous system of 1st order ODEs
%   
%   Compare numerical solution using ODE23 with analytical solution.
%   See ODE function file Tanks2_unforced_eqns.m for equations evaluated.
%
%   In this case we have a sequence of two interconnected tanks.  There
%   is no salt entering the system (i.e. unforced system).  Tank 1 has some
%   initial salt content but the other one does not.  Tank 1 feeds Tank 2 with
%   a constant flow rate in the system (i.e. tank volumes remain constant).
%   The goal is to determine the amount of salt in each tank versus time.  
%   The variables of interest are:
%     m(i)  = mass of salt in tank i (lbm) (numerical)
%     me(i) = mass of salt in tank i (lbm) (exact)
%     v(i)  = volume of tank i (gal)
%
%   File prepared by J. R. White, UMass-Lowell (Sept. 2003)
%

%
%   getting started
      clear all,  close all
%
%   Define basic variables
      global v q 
      v = [20 40];         % fixed volume of tanks  (gal)
      mo = [15 0];         % initial mass of salt in tanks (lbm)
      q = 10;              % constant flow rate in system (gal/min)
%
%   Now solve 2x2 system of ODEs (Initial Value Prob -- IVP)
      to = 0;     tf = 25;
      [t,m] = ode23('tanks2_unforced_eqns',[to tf],mo);
%
%   Exact soln to sample problem (uses specific data from hand calc - not general)
      te = linspace(to,tf,51)';   me = zeros(length(te),2);
      lam = [-0.5 -0.25];  x1 = [1 -2]';   x2 = [0 1]';  
      c1 = 15;   c2 = 30;   
      me(:,1) = c1*x1(1)*exp(lam(1)*te);
      me(:,2) = c1*x1(2)*exp(lam(1)*te) + c2*x2(2)*exp(lam(2)*te);
%
%   Compare Exact and Numerical Solution for sample problem
      figure(1)
      plot(te,me(:,1),'g-',te,me(:,2),'g--', ...
           t,m(:,1),'ro',t,m(:,2),'r+','LineWidth',2)
      title('Tanks2\_Unforced:  Two Mixing Tanks in Series (Unforced)')
      xlabel('time (min)'),ylabel('mass of salt (lbm)'),grid
      gtext('tank 1'),gtext('tank 2')
      ss = ['Exact -- lines     ';'Numerical -- points']; ss = {ss};
      gtext(ss)
%
%   End of Problem
%
