%% Math 4480/5480 - Example 3.1. Nonlinear Lamb (p.61-71)

%%
clear all; close all;
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*exp(0.025*x)-0.45*x;
ezplot(y,[0,40]);

%%
% Differentiate the symbolic expressions

dydx = diff(y,x)

%%
% Find the critical points of y=y(x). The *solve* command attempts to find explicit expressions for
% the roots of the equation dydx=0. This is in many cases impossible. 
% If symbolic solutions cannot be determined, numeric solutions are
% returned. (See below)

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

%% 
% Notice that we get two critical points (roots of dydx=0). We must pick only values that have a
% meaning in our original problem, in this case the positive value only.

xmax =xmax(1)

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

ymax=subs(y,x,xmax)

%%
% Alternately, numeric solutions of dydx=0 can be obtained through *Newton's Method*: 

F = dydx;
F1 = diff(F,x);

format long
N = 10; % number of iterations
x0 = 19 % initial guess
fprintf(' iteration  xvalue\n\n');
for i=1:N
    x1=x0-subs(F,x,x0)/subs(F1,x,x0);
    fprintf('%5.0f       %1.16f\n', i, x1);
    x0 = x1;
end

%%
display('Hence, the critical point (solution of F=0) is (approx)'), x1


%%
% Sensitivity to the parameter *r* (p. 10)

%%
% It turns out that, in this particular case, one could AGAIN solve 'explicitly' 
% the equation dydx=0 when *r* is treated as a symbolic variable. 
% This is though plain luck, since there exists a special function exists (and is
% programmed in the computer) that describes the solution of dydx=0.
% Attempting to go explicit when studying sensitivity is highly NOT RECOMMENDED way to go!

%%
% Sensitivity to the parameter *r* (p. 10)
 xvalues = 0;
 for r = 0.008:0.001:0.012
     y = (0.65-r*x)*200*exp(0.025*x)-0.45*x;
     dydx=diff(y,x);
     xmaxr=solve(dydx,x);
      xmaxr = double(xmaxr);
      xmaxr = xmaxr(1);
     xvalues = [xvalues; xmaxr];
 end
 xvalues = xvalues(2:end);
 
%%
% Display the optimal times x for several values of c around 0.025 
 rvalues = [0.008; 0.009; 0.01; 0.011; 0.012];
 format short; 
 display([rvalues,xvalues])

%% 
% Plot the xmax as a function of r, indicating the values shown in the table above 
 figure, plot(rvalues,xvalues, '--mo',...
                'LineWidth',2,...
                'MarkerEdgeColor','k',...
                'MarkerFaceColor',[.49 1 .63],...
                'MarkerSize',12)
xlabel('r'), ylabel('xmax')


%%
% Sensitivity to the parameter *c* (p. 10)
 xvalues = 0;
 for c = 0.022:0.001:0.028
     y = (0.65-0.01*x)*200*exp(c*x)-0.45*x;
     dydx=diff(y,x);
     xmaxc=solve(dydx);
      xmaxc = double(xmaxc);
      xmaxc = xmaxc(1);
     xvalues = [xvalues; xmaxc];
 end
 xvalues = xvalues(2:end);
 
%%
% Display the optimal times x for several values of c around 0.025 
 cvalues = 0.022:0.001:0.028;
 cvalues=cvalues'; % transposes the row into a column
 format short; display([cvalues,xvalues])
 

%% 
% Plot the xmax as a function of r, indicating the values shown in the table above 
 figure, plot(cvalues,xvalues, '--mo',...
                'LineWidth',2,...
                'MarkerEdgeColor','k',...
                'MarkerFaceColor',[.49 1 .63],...
                'MarkerSize',12)
xlabel('c'), ylabel('xmax')
