function x = multiplefsolve(perimeter,area,jacflag);
%x = multiplefsolve(perimeter,area,jacflag);
%find the sides of a rectangle with given perimeter and area
%note that there may be NO SOLUTION in cases where the one fourth of the 
%perimeter is less than the square root of the area...why?
%
%VARIABLES:
%perimeter: perimeter of desired rectangle
%area: area of desired rectangle
%jacflag: set to 1 if want to calculate jacobian matrix explicitly, else 0
%
%EXAMPLE:
%consider a 4 x 3 rectangle...
%perimeter = 14, area = 12, jacflag = 1 or 0;
%x = multiplefsolve(14,12,1); 
%x = multiplefsolve(14,12,0);

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


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


return;



function [y,Jacobian] = g_with_jacobian(x,param);
    L = x(1); %rectangle length
    W = x(2); %rectangle width
    
    y1 = 2*L + 2*W - param.perimeter;
    y2 = L*W - param.area;
    
    dy1dL = 2;
    dy1dW = 2;
    dy2dL = W;
    dy2dW = L;
    
    Jacobian = [ dy1dL dy1dW ;...
                 dy2dL dy2dW ];
    
    y = [y1 y2]';
return;

%if jacobian is set to off, then then the function is just:
function y = g_no_jacobian(x,param);
    L = x(1); %rectangle length
    W = x(2); %rectangle width
    
    y1 = 2*L + 2*W - param.perimeter;
    y2 = L*W - param.area;
    
    y = [y1 y2]';
return;