function [R,H] = sampleproblem5(area,volume,jacflag,R0,H0);
%[R,H] = sampleproblem5(area,volume,jacflag,R0,H0);
%find radius and height of a cylinder with given surface area and volume
%note that there may be NO SOLUTION or MULTIPLE SOLUTIONS
%
%VARIABLES:
%area: surface area of desired cylinder
%volume: volume of desired cylinder
%jacflag: set to 1 if want to calculate jacobian matrix explicitly, else 0
%R0,H0: initial guesses for radius and height
%EXAMPLE:
%consider a cylinder with R = 1 and H = 1...
%area = 2*(pi*R^2) + (2*pi*R*H) = 4pi, volume = pi*R^2*H = pi;
%SOUGHT SOLUTION
%[R,H] = sampleproblem5(4*pi,pi,0,1.1,0.9)
%[R,H] = sampleproblem5(4*pi,pi,1,1.1,0.9)
%ANOTHER SOLUTION
%[R,H] = sampleproblem5(4*pi,pi,1,0.5,3)
%Initial guesses are important!

%BEGIN CODE
x0 = [R0,H0];%[volume^1/3,sqrt(area)]'; %initial guess
tol = sqrt(eps)/100; %sqrt(eps) ~ 10^-8, often a reasonable tolerance
param.volume = volume;
param.area = area;
format long;


if(jacflag == 0)    
    options = optimset('TolFun',tol,'TolX',tol,'Jacobian','off');
    x = fsolve(@gcyl,x0,options,param);
else
    options = optimset('TolFun',tol,'TolX',tol,'Jacobian','on');
    x = fsolve(@gcyl_jac,x0,options,param);
end

R = x(1);
H = x(2);
return;



function [y,Jacobian] = gcyl_jac(x,param);
    R = x(1); %cylinder radius
    H = x(2); %cylinder height
    
    y1 = 2*pi*R^2 + 2*pi*R*H - param.area;
    y2 = pi*R^2*H - param.volume;
    
    dy1dR = 4*pi*R + 2*pi*H;
    dy1dH = 2*pi*R;
    dy2dR = 2*pi*R*H;
    dy2dH = pi*R^2;
    
    Jacobian = [ dy1dR dy1dH ;...
                 dy2dR dy2dH ];
    
    y = [y1 y2]';
return;

%if jacobian is set to off, then then the function is just:
function y = gcyl(x,param);
    R = x(1); %cylinder radius
    H = x(2); %cylinder height
    
    y1 = 2*pi*R^2 + 2*pi*R*H - param.area;
    y2 = pi*R^2*H - param.volume;

    y = [y1 y2]';
return;