How to combine nested functions?

I am writing a program where some functions are nested into a parent function. I have a similar code in gauss which gets the job done easily but I wish to understand how to do a similar exercise in matlab.

Essentially, I have a parent function which performs a golden section search to a sub-function which is, in my case, a value function (a bellman equation) defined over a domain x.

In the base environment in gauss I call:

aropt[l,i]= golden(&value1,ax,bx,cx,tol1);

where “golden”, the parent function, is a standard golden search algorithm given by:

proc golden(&f,ay,by,cy,tol);
    local f:proc,x0,x1,x2,x3,xmin,r1,r2,f1,f2;
    r1=0.61803399;
    r2=1-r1;
    x0=ay;
    x3=cy;
    if abs(cy-by)<=abs(by-ay);
        x1=by;
        x2=by+r2*(cy-by);
    else;
        x2=by;
        x1=by-r2*(by-ay);
    endif;
    f1=-f(x1);
    f2=-f(x2);
    do until abs(x3-x0)<=tol*(abs(x1)+abs(x2));
        if f2<f1;
            x0=x1;
            x1=x2;
            x2=r1*x1+r2*x3;
            f1=f2;
            f2=-f(x2);
        else;
            x3=x2;
            x2=x1;
            x1=r1*x2+r2*x0;
            f2=f1;
            f1=-f(x1);
        endif;
    endo;
    if f1<=f2;
        xmin=x1;
    else;
        xmin=x2;
    endif;
    retp(xmin);
endp;

And that handles a nested function f (i.e. “value1”) defined over x, containing both input variables and parameters and variables defined internally.

proc value1(x);
    local c;
    c=(1+r)*a[l]+pen-x;
    if c<=0;
        retp(-1e10);
    endif;
    retp(u(c,1)+b*rvalue(x,i+1));
endp;

Which in turn read other sub functions such as “u” and “rvalue”.

To replicate this into matlab I followed the basic documentation on nested functions( https://www.mathworks.com/help/matlab/matlab_prog/nested-functions.html) but so far I have been not successful in doing so.

As far I understood, in I have to pass in the whole bunch of variables and parameters used as inputs in nested functions also to their parent functions, unless I specify them as global variables.

Hence, this would leave me with

p = GSVF1(aa,bb);
aropt(l,i)=p(1);

which returns:

Not enough input arguments.
Error in value1 (line 3)
cr= (1+r)*a(l)+pen-a(m);
Error in goldensearchmax (line 24)
fl = f(xl);
Error in HM91v_A18small (line 379)
[z] = goldensearchmax(@value1,ax,cx);

In the code above a(m) is the input variable for the value function (corresponding to x in “value1” in gauss). To simplify things, I incorporated the utility and residual function functions directly in the body of value1 but I am not sure that I can do it.

function [p] = GSVF1(aa,bb)
% -------------------------------------------------------------------------
% Draft program that incorporates a maximization problem into a golden
% section search algorithm. More specifically, the program tests how to
% nest a value function and the dependent optimal policy function for
% capital into a parent function that performs a golden search algorithm.
% Valid for retired agents only.
%--------------------------------------------------------------------------
% INPUTS
%   (f: Declared function in base environment.)
%   a: Lower bound of the initial interval to be evaluated
%   b: Upper bound of the initial interval to be evaluated
%
% OUTPUTS
%   m = [m1 m2], where m1 is the maximum and m2 is the maximizer
%--------------------------------------------------------------------------
p = @value1;
    function [VFr] = value1(a,m,r,l,pen,beta,i,kmax,na,vr,VI_method,sigma,gamma,psi,empl)
        %------------------------------------------------------------------
        % INPUTS
        % - parameters:             beta, sigma, gamma, psi, pen
        % - factor prices:          r
        % - asset grid:             a, na, kmax, l, m
        % - age:                    i
        % - labor supply (=0):      empl
        % - other:                  VI_method
        %
        % OUTPUTS
        %   p = [VFr] Value function for the retired agent at age i with asset l.
        %------------------------------------------------------------------
        cr= (1+r)*a(l)+pen-a(m);
        if cr <= 0
            VFr   = -1e10;
        elseif cr > 0
            % Utility -----------------------------------------------------
            if sigma == 1
                utilr = log(cr+psi)+gamma*log(1-empl);
            else
                utilr = ((((cr+psi)*(1-empl)^gamma)^(1-sigma)-1))/(1-sigma);
            end
            % Residual Value ----------------------------------------------
            t0 = i+1;
            k  = a(m);
            k0 = k/kmax*(na-1)+1;
            n2 = floor(k0);
            n1 = k0-n2;
            
            if k <= 0
                RVr = vr(1,t0)-(1-k0)*(vr(2,t0)-vr(1,t0));
            elseif k >= kmax
                RVr = vr(na,t0);
            else
                if VI_method == 1 % linear interp
                    RVr = (1-n1)*vr(n2,t0)+n1*vr(n2+1,t0);
                elseif VI_method == 2 % cublic spline
                    disp ('error: method 2 not allowed - change to method 1')
                end
            end
            % Value Function ----------------------------------------------
            VFr = utilr+beta*RVr;
        end
    end
 
% Golden Search -----------------------------------------------------------
xl   = aa+(3-5^.5)/2*(bb-aa);
xh   = aa+(5^.5-1)/2*(bb-aa);
fl   = p(xl);
fh   = p(xh);
d    = abs(xh-xl);
tol  = 10^-10;
iter = 0;
 
while d > tol
    iter=iter+1;
    d = abs(xl-xh);
    if f(xl)<f(xh)
        aa = xl;
        xl = aa+(3-5^.5)/2*(bb-aa); fl = f(xl);
        xh = aa+(5^.5-1)/2*(bb-aa); fh = f(xh);
    else
        bb = xh;
        xl = aa+(3-5^.5)/2*(bb-aa); fl = f(xl);
        xh = aa+(5^.5-1)/2*(bb-aa); fh = f(xh);
    end
end
 
if fh>fl
    p = [xh];
else
    p = [xl];
    
end
end

I would appreciate any help on how to combine parent and nested functions with function handles. I believe there are simpler ways to achieve this but to start with any hint would be great.

Thank you,

Andre