%   
% Written by:
% -- 
% John L. Weatherwax                2006-08-28
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

close all; clc; clear; 

% set the constants of this problem: 
t_0     = 0.; 
t_f     = 12.; 
xair    = 45;
xd      = 70; 
x0      = 45; % the initial conditions on x_1(t) 
T       = 4;
k1      = 1.0; 
c1      = 5.0; c2 = 1.0; % the ratio needs to be 5.

% to get the initial guess for the fzero function we plot the right and left hand sides of the equation we want to solve and look where they cross: 
% 
x1f = linspace( 45, 70, 100 );
rhs_function = @( x1f ) ( xair + (x0 - xair) * exp( -(t_0+t_f)/T ) - ( (k1^2)*T*(x1f - xd)*c1/(2*c2) ) * ( 1 + exp( 2*(t_0-t_f)/T ) ) );

rhsh = plot( x1f, rhs_function( x1f ), '-bx' ); hold on;  
lhsh = plot( x1f, x1f, '-r' );
xlabel('x1f'); legend( [rhsh,lhsh], {'right-hand-side','left-hand-side'} );

x1f_guess = 67.73; % = a value that is very close to the root
x1f_zero = fzero( @( x )( rhs_function(x) - x ), x1f_guess )

% we can now plot x_1(t) and u(t) using the known analytic representations of them as desired
%
t = linspace( t_0, t_f, 100 ); 

x1 = xair + ( x0 - xair ) * exp( -(t + t_0)/T ) - ( (k1^2)*T*(x1f_zero - xd)*c1/(2*c2) ) * ( exp( (t - t_f)/T ) + exp( 2*(t_0-t_f)/T ) );
figure; plot( t, x1, '-o' ); xlabel('time'); ylabel('x_1(t)'); 
%saveas( gcf, '../../WriteUp/Graphics/Chapter3/sect_5_prob_4_x1.eps', 'epsc' ); 

u = - (k1*c1/c2)*(x1f_zero - xd)*exp( (t-t_f)/T );
figure; plot( t, u, '-o' ); xlabel('time'); ylabel('u(t)'); 
%saveas( gcf, '../../WriteUp/Graphics/Chapter3/sect_5_prob_4_u.eps', 'epsc' );