function SolveContact

global v d h eo po n

% Define known inputs
v = 0.45;
d = 0.1;
h = 8;
eo = 0.06;
po = 18.4;
n = 3;

% Plot E as a function of p to estimate solution
E = linspace(0,1000,101)';
p = linspace(0,50,101)';

pcontact = ContactEquation(E);
Ematerial = MaterialEquation(p);

plot(pcontact,E,'k-',p,Ematerial,'r-')
xlabel('p')
ylabel('E')
pause

% Approach 1: Substitute material equations into pressure equation
% and solve one nonlinear equation for p. Then substitute p value into
% material equation and calculate for E.

fprintf('Using rootfinding to solve one nonlinear equation in one unknown . . .\n')
p0 = 20;
options = optimset('Display','iter');
p = fzero(@OneZeroEqn,p0,options)
E = MaterialEquation(p)
pause

% Approach 2: Solve two nonlinear equations simultaneously for p and E.
fprintf('Using rootfinding to solve two nonlinear equations in two unknowns . . .\n')
x0(1,1) = 20;
x0(2,1) = 300;
options = optimset('Display','iter');
x = fsolve(@TwoZeroEqns,x0,options);
p = x(1,1)
E = x(2,1)
pause

% Approach 3: Solve unconstrained nonlinear optimization problem for p and E
% Perform the optimization using three different unconstrained optimization
% algorithms
fprintf('Using nonlinear optimization to solve for two design variables . . .\n')
x0(1,1) = 20;
x0(2,1) = 300;
lb = [];
ub = [];
options = optimset('Display','iter','MaxFunEvals',500);
fprintf('First try lsqnonlin optimization algorithm . . .\n')
x = lsqnonlin(@CostFunc1,x0,lb,ub,options); % This one does not work
p = x(1,1)
E = x(2,1)
pause
fprintf('Next try fminsearch optimization algorithm . . .\n')
x = fminsearch(@CostFunc2,x0,options); % This one works
p = x(1,1)
E = x(2,1)
pause
fprintf('Finally try fminunc optimization algorithm . . .\n')
x = fminunc(@CostFunc2,x0,options); % This one works
p = x(1,1)
E = x(2,1)

%-------------------------------------------------------------------------
function p = ContactEquation(E)

global v d h eo po n

p = ((1-v)/((1+v)*(1-2*v)))*E*(d/h);

%-------------------------------------------------------------------------
function E = MaterialEquation(p)

global v d h eo po n

E = 1./((0.5*eo/po)*(1+n*(p/po).^(n-1)));

%-------------------------------------------------------------------------
function zero = OneZeroEqn(p)

global v d h eo po n

E = 1./((0.5*eo/po)*(1+n*(p/po).^(n-1)));
zero = ((1-v)/((1+v)*(1-2*v)))*E*(d/h)-p;

%-------------------------------------------------------------------------
function zero = TwoZeroEqns(x)

global v d h eo po n

p = x(1,1);
E = x(2,1);

zero = zeros(2,1);
zero(1,1) = ((1-v)/((1+v)*(1-2*v)))*E*(d/h)-p;
zero(2,1) = 1/((0.5*eo/po)*(1+n*(p/po)^(n-1)))-E;

%-------------------------------------------------------------------------
function errors = CostFunc1(x)

global v d h eo po n

p = x(1,1);
E = x(2,1);

errors = zeros(2,1);
errors(1,1) = ((1-v)/((1+v)*(1-2*v)))*E*(d/h)-p;
errors(2,1) = 1/((0.5*eo/po)*(1+n*(p/po)^(n-1)))-E;

%-------------------------------------------------------------------------
function error = CostFunc2(x)

errors = CostFunc1(x);

error = errors(1,1)^2+errors(2,1)^2;