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

clc; 
close all; 

% solve for the value of x_P:
rt_fn = @(x) (-0.5*pi*(1+x).*sin(pi*x) - 0.5 + sin(0.5*pi*x).^2); 
x_P = fsolve( rt_fn, -0.1 ); 

% solve for the value of x_Q: 
rt_fn = @(x) (-0.5*pi*(-1+x).*sin(pi*x) - 0.5 + sin(0.5*pi*x).^2); 
x_Q = fsolve( rt_fn, +0.1 ); 

all_x = linspace(-1,+1);

% get the coefficients of the linear function to the left: 
B = -0.5*pi*sin(pi*x_P); 
A = B; 

% get the coefficients of the linear function to the right:
D = -0.5*pi*sin(pi*x_Q); 
C = -D; 

indL = find( all_x <= x_P ); 
indM = find( (all_x > x_P) & (all_x < x_Q) ); 
indR = find( all_x >= x_Q ); 

% map each region into its functional representation: 
obs_y = zeros(1,length(all_x)); 
obs_y(indL) = A+B*all_x(indL); 
obs_y(indM) = 0.5 - sin(0.5*pi*all_x(indM)).^2; 
obs_y(indR) = C+D*all_x(indR); 

figure; sh=plot( all_x, obs_y, '-k' ); hold on; grid on; 
        oh=plot( all_x, 0.5-sin(0.5*pi*all_x).^2, '-go' ); 
axis( [-1,+1,0,Inf] ); legend([sh,oh],{'obstacle solution','obstacle'},'location','south'); 
xlabel('x-axis'); ylabel('u(x)');
title( 'obstacle solution' );
saveas( gcf, '../../WriteUp/Graphics/Chapter7/prob_7_3_pt_b', 'epsc' );