%
%   Tanks2_forced1.m    Numerical Integration of System of ODEs (IVP)
%
%   Sample problem from in-class notes -- forced system of 1st order ODEs
%   
%   Compare numerical solution using ODE23 with analytical solution.
%   See ODE function file Tanks2_forced1_eqns.m for equations evaluated.
%
%   In this case we have a sequence of two interconnected tanks.  Salt enters
%   the system (i.e. forced system) with Ci = 0.5exp(-0.25t).  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 (original Oct. 2001)
%

%
%   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_forced1_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 = -1;   c2 = 3;   
      me(:,1) = 5*c1*exp(lam(1)*te) + 20*exp(lam(2)*te);
      me(:,2) = -10*c1*exp(lam(1)*te) + (10*te - 40 + 10*c2).*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\_Forced1:  Two Mixing Tanks in Series (Forced -- g_1(t) = 5e^{-0.25t})')
      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
