Returning multiple values in fminbnd

I am trying to solve a dynamic optimization problem in which I solve for optimxal capital and debt. The problem is structured such that for each period I have to find the value of debt, variable d, that maximizes the firm value, given the value of k and k0. And then in the next step I have to find the value of k that maximizes the firm value over time, given k0. My main looks as follows:

disc  = 0.96; %discount rate for the expected next period value function
alpha = 0.45; %capital elasticity, Cobb-Douglas *
beta  = 0.34; %labor elasticity, Cobb-Douglas *

R     = 719.0615; %Market Size
f     = 4.02; %fixed cost of production
sigma = 5.23; %elasticity of substitution
W     = 0.3;  %wage *
theta = 0.7;  %cost of bankruptcy *
chi   = 0.8;  %level of debt at which bankruptcy happens *
c     = 0.2;  %parameters of capital adjustment *
kappa = 0.3;  %managerial rent *
mu_b  = 0.1;  %mean of the firm level shock * 
sigma_b = 0.12; %standard deviation of the firm level shock 

omega = (beta*(sigma-1))/(beta*(sigma-1) - sigma);%an itermediary parameter

na = 3;%TFP shock grid size
nk = 11; %capital grid size

%setting up the capital grid
kmin = 1;
kmax = 10000000000;
kvec = linspace(log(kmin),log(kmax),nk)';
kvec = exp(kvec);

%setting up the TFP grid
multiple = 3;
mu_a     = 0;
rho_a    = 0.7;
sigma_a  = 0.05;
[avec,pi]=tauchen(0.0,sigma_a,rho_a,multiple,na);

Unrecognized function or variable 'tauchen'.

avec = exp(avec)';


tol = 1e-5;
diff = 1;
its  = 0;
maxits = 500;

v0 = zeros(nk,na);

while diff>tol && its <maxits
    for j = 1:na
        for i = 1:nk
            P0 = 0.001;
            a0 = avec(j);
            k0 = kvec(i);
            k1 = fminbnd(@(k)ValueFunctionCapital(k,v0,k0,a0,j,P0),kmin,kmax);
            [val,d] = -ValueFunctionCapital(k1,v0,k0,a0,j,P0);
            P = fixedPointP(P0,k1);
            
            v1(i,j) = val;
            k11(i,j) = k1;
            d11(i,j) = d;
            p11(i,j) = P;
            
            
        end 
    end
    diff = norm(v0-v1)
    v0 = v1;
    its = its+1
end

My ValueFuntionFunctionCapital looks as follows:

function [val,d] = ValueFunctionCapital(k,v0,k0,a0,j,P)
 
global disc alpha beta R f sigma W theta chi c kappa mu_b sigma_b omega pi mu_a sigma_a kvec avec maxits
 
   g1 = interp1(kvec,v0,k,'linear'); %smooth out previous value function
   
   %some constants 
   prof_brac = (sigma/(beta*(sigma-1)))^omega - (sigma/(beta*(sigma-1)))...
         ^(sigma/(beta*(sigma-1)-sigma));
   quad_cost = c*(k0-k)^2;
   function val_inner = ValueFunctionDebtInner(d,k0,k,a0,P)
       g_inner = @(b) lognpdf(b, mu_b, sigma_b); % define PDF of lognormal distribution
       fun_inner = @(b) (1./b).^omega .* g_inner(b);
       Q11_inner = integral(fun_inner, 0, Inf);
 
       pdf_inner = @(x) lognpdf(x, mu_b, sigma_b);
 
       V_term1_1_inner = (1-kappa)*(R*P^(sigma-1))^(-1/(beta*(sigma-1)))*(a0*k^alpha)^...
            (-(sigma-1)/(beta*(sigma-1)-sigma))*W^omega*prof_brac*Q11_inner;
       V_term1_2_inner = (1-kappa)*(k0-k-f*k-quad_cost);
       V_term1_inner = V_term1_1_inner + V_term1_2_inner;
   
   
       b_th_inner = bankruptcy(a0,k0,k,d,P,prof_brac,W,chi,kappa,f,sigma,beta,c,R);
       V_term2_cdf_inner = logncdf(b_th_inner, mu_b, sigma_b); % compute CDF at x ;
       V_term2_inner = kappa*d*(1-V_term2_cdf_inner);
   
 
       Q31_inner = integral(fun_inner, 0, b_th_inner);
   
      % define PDF function using lognpdf
       Q32_inner = integral(pdf_inner, 0, b_th_inner);
   
       V_term3_1_inner = (1-kappa)*theta*(R*P^(sigma-1))^(-1/(beta*(sigma-1)-sigma))*...
           (a0*k^alpha)^ (-(sigma-1)/(beta*(sigma-1)-sigma))*W^omega*...
           prof_brac*Q31_inner;
       V_term3_2_inner = (1-kappa)*theta*(k0-k-f*k-quad_cost)*Q32_inner;     
       V_term3_inner = V_term3_1_inner + V_term3_2_inner;
    
 
       V_final_inner = V_term1_inner + V_term2_inner - V_term3_inner ;
          
       val_inner = V_final_inner;
       val_inner = -val_inner;
       
       
   end
 
   dmin = 1;
   dmax = k;
   
   d = fminbnd(@(d)ValueFunctionDebtInner(d,k0,k,a0,P),dmin,dmax);
 
 
   g = @(b) lognpdf(b, mu_b, sigma_b); % define PDF of lognormal distribution
   fun = @(b) (1./b).^omega .* g(b);
   Q11 = integral(fun, 0, Inf);
 
   pdf = @(x) lognpdf(x, mu_b, sigma_b);
 
   V_term1_1 = (1-kappa)*(R*P^(sigma-1))^(-1/(beta*(sigma-1)))*(a0*k^alpha)^...
            (-(sigma-1)/(beta*(sigma-1)-sigma))*W^omega*prof_brac*Q11;
   V_term1_2 = (1-kappa)*(k0-k-f*k-quad_cost);
   V_term1 = V_term1_1 + V_term1_2;
   
   
   b_th = bankruptcy(a0,k0,k,d,P,prof_brac,W,chi,kappa,f,sigma,beta,c,R);
   V_term2_cdf = logncdf(b_th, mu_b, sigma_b); % compute CDF at x ;
   V_term2 = kappa*d*(1-V_term2_cdf);
   
 
   Q31 = integral(fun, 0, b_th);
   
    % define PDF function using lognpdf
   Q32 = integral(pdf, 0, b_th);
   
   V_term3_1 = (1-kappa)*theta*(R*P^(sigma-1))^(-1/(beta*(sigma-1)-sigma))*...
           (a0*k^alpha)^ (-(sigma-1)/(beta*(sigma-1)-sigma))*W^omega*...
           prof_brac*Q31;
    V_term3_2 = (1-kappa)*theta*(k0-k-f*k-quad_cost)*Q32;     
    V_term3 = V_term3_1 + V_term3_2;
    %V_term3 = -V_term3;
 
    V_final = V_term1 + V_term2 - V_term3 ;
          
    val = V_final + disc * (g1*pi(j,:)');
    val = -val;
end

I understand that there is a lot of repetition of code here. But for now what I want is that I can store k, the value of capital that maximizes the value and the corresponding d. I am having troubles storing that since with fminbnd I am not able to retrieve both the values when I call the function. Is there anyway I can get both k and corresponding d in main?

[Torsten]

It's not clear to me at which position in your code you want to output what quantities.

[Prerna Mishra]

k1 = fminbnd(@(k)ValueFunctionCapital(k,v0,k0,a0,j,P0),kmin,kmax);

In the above line I want both the optimized k and d.

[Torsten]

You can get d only by calling "ValueFunctionCapital" again with the optimal k1 as you do in the next line of your code. Doesn't that work ? I cannot run your code since at least "tauchen" is missing.