%% MATH 448/548 - Example 1.1. The Lamb Problem (p.3-14)
% This is a Matlab code which illustrates some of the symbolic capabilities
% of Matlab reflected in Example 1.1 of your textbook. 
% The HTLM version is generated with the useful publish command.

%% Introduce the symbolic (independent) variable(s)
% It is generally recommended to *clear all* previously defined
% variables from the workspace. Use *clf* to clear all previously defined
% figures.

clear all; clf;
syms x

%% Define (and plot) the function (dependent variable)
% Symbolic functions (of one or several symbolic variables) 
% can be defined and plotted as follows. 

y = (0.65-0.01*x)*(200+5*x)-0.45*x;
ezplot(y,[0,20]);

%% Differentiate symbolic expressions
% Matlab can perform symbolic differentiation using *diff* command

dydx = diff(y,x);

%% Find the critical points of y=y(x), 
% This is done by solving the equation dydx=0. So far xmax is still thought 
% as symbolic variable. Matlab can convert symbolic into 
% numeric (double) format

xmax = solve(dydx);
xmax = double(xmax)

%% Compute the numerical value ymax=y(xmax)

ymax=subs(y,x,xmax);

%% Plot the symbolic function, with specified domain 
% Use *ezplot* command to plot symbolic functions (notice the syntax for specifying  
% the domain of the function, in this case [0,20]). Use *plot* command for numeric
% functions or variables. Commands *hold on* and *off* allow for plotting
% different plots on the same figure. 

ezplot(y,[0 20])
hold on
plot(xmax,ymax,'-ro');
grid on
hold off


%% Sensitivity to the parameter *r* (p. 10)
% Introduce variable price drop *r* and compute the sensitivity

syms r
y = (0.65-r*x)*(200+5*x)-0.45*x;
dydx=diff(y,x);
xmaxr=solve(dydx,x)

%% Display optimal times to sell for several values of *r* (see Table 1.1
%% on page 10)
rvalues=[0.008; 0.009; 0.01; 0.011; 0.012];
xvalues = subs(xmaxr, r, rvalues);
display([rvalues, xvalues])

%% Plot the xmax as a function of r, indicating the values shown in the table above 
figure, ezplot(xmaxr,[0.008 0.012])
hold on
plot(rvalues,xvalues, 'ks')

%%Compute the sensitivity of x to r, S(x,r) when r=0.01
Sr = diff(xmaxr,r)*r/xmaxr;
S = subs(Sr, r, 0.01) 

%% Sensitivity to the parameter *g* (p. 13)
% Introduce variable growth rate g 

syms g
y = (0.65-0.01*x)*(200+g*x)-0.45*x;
dydx=diff(y,x);
xmaxg=solve(dydx,x);

%% Plot of the optimal time vs g
% The command *figure* instructs Matlab to open a new figure for potting

figure, ezplot(xmaxg,[3,7]);

%% Evaluate sensitivity 
% This is the sensitivity of xmaxg to variations of g aroung the value g=5

dxdg=diff(xmaxg,g);
Sxg=dxdg*g/xmaxg;
S=subs(Sxg,g,5)
