function Show_fsolve
% Illustrates several canonical applications of fsolve which
% can be used to find zeros of f:Rn--->Rn.
% See the functions MyF1, MyF2, and MyF3 below. They are used
% in the examples.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Example 1.

% Finds a zero of MyF1. 
% No Jacobian Supplied.
% Uses default tolerances.
% Returns just the zero.

clc
disp('Example 1.')
tInitial = randn(2,1);
% Turns off messaging...
options = optimset('Display','off');
t_star  = fsolve(@MyF1,tInitial,options)
pause

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Example 2.

% Finds a zero of MyF2. 
% Jacobian Supplied.
% Uses default tolerances.
% Returns just the zero.

clc
disp('Example 2.')
tInitial = randn(2,1);
options = optimset('Display','off','Jacobian','on');
t_star  = fsolve(@MyF2,tInitial,options)
pause

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Example 3.

% Finds a zero of MyF1. 
% No Jacobian Supplied.
% Uses default tolerances.
% Returns a zero and the value of F at the zero

clc
disp('Example 3.')
tInitial = randn(2,1);
options = optimset('fsolve'); % reset to default settings
options = optimset('Display','off');
[t_star,MyF1_star] = fsolve(@MyF1,tInitial,options)
pause

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Example 4.
% Finds a zero of MyF3(t,E1,E2). 
% No Jacobian Supplied.
% Uses default tolerances.
% Returns a zero and the value of MyF3 at the zero
% Shows how to apply fsolve to a function with parameters

clc
disp('Example 4.')
tInitial = randn(2,1);
E1 = struct('P',3,'A',10,'tau',pi/6);
E2 = struct('P',4,'A',7,'tau',-pi/3);
options = optimset('Display','off');
[t_star,MyF3_star] = fsolve(@(t) MyF3(t,E1,E2),tInitial,options)
pause

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Example 5.
% Finds a zero of MyF1(t). 
% No Jacobian Supplied.
% Sets some useful tolerances.
% Returns a zero and the value of MyF1 at the zero

clc
disp('Example 5.')
tInitial = randn(2,1);

% fsolve will quit as soon as any of these criteria are met...
del_x = .00000000001; 
del_y = .001;  
fEvalMax = 100;
itMax = 20;
options = optimset('MaxFunEvals',fEvalMax,'MaxIter',itMax);
options = optimset('TolX',del_x,'TolFun',del_y);
options = optimset('Display','off');

[t_star,MyF1_star] = fsolve(@MyF1,tInitial,options)
pause

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Example 6.
% Finds a zero of MyF3(t,E1,E2). 
% No Jacobian Supplied.
% Sets tolerances.
% Returns a zero, the value of MyF3 at the zero, and output information.
% Shows how to apply fsolve to a function with parameters

clc
disp('Example 6.')
tInitial = randn(2,1);
E1 = struct('P',3,'A',10,'tau',pi/6);
E2 = struct('P',4,'A',7,'tau',-pi/3);
[t_star,MyF1_star,ExitFlag,Output] = fsolve(@(t) MyF3(t,E1,E2),tInitial,options);

TheComputedRoot =  t_star
FunctionAtRoot  =  MyF1_star
fprintf('Iterations  = %1d\n',     Output.iterations)
fprintf('Fevals      = %1d\n',     Output.funcCount)
fprintf('Method      = %s\n',      Output.algorithm)
if ExitFlag==1
   fprintf('Exit Status = FSOLVE converged to a solution\n')
elseif ExitFlag==2
   fprintf('Exit Status = Change in X <  tol\n')
elseif ExitFlag==3
   fprintf('Exit Status = Change in the residual < tol\n')
elseif ExitFlag==4
   fprintf('Exit Status = Search direction < tol\n')
elseif ExitFlag ==0
   fprintf('Exit Status = Too many iterations/fEvals\n')
elseif ExitFlag == -2
   fprintf('Exit Status = Converging to non-root.\n')
elseif ExitFlag == -3
   fprintf('Exit Status = Trust region radius became too small.\n')
else
   fprintf('Exit Status = Line search not good enough.\n')
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Functions used in the examples....

function y = MyF1(t)
E1 =  [5*cos(t(1)); 3*sin(t(1))];
E2 = [cos(pi/6) sin(pi/6); -sin(pi/6) cos(pi/6)]*[cos(t(2)); 10*sin(t(2))];
y = E1-E2;

function [y,J] = MyF2(t)
y = MyF1(t);
% Compute a finite difference Jacobian
e1 = [1;0]; e2 = [0;1]; h = .00000001;
J = [(MyF1(t+h*e1)-y)/h (MyF1(t+h*e2)-y)/h];


function z = MyF3(t,E1,E2)
z = Orbit(t(1),E1)-Orbit(t(2),E2);

function z = Orbit(t,E)
Q = [cos(E.tau) -sin(E.tau); sin(E.tau) cos(E.tau)];
x = (E.P-E.A)/2 + ((E.P+E.A)/2)*cos(t);
y = sqrt(E.P*E.A)*sin(t);
z = Q*[x;y];


 