c_lo = 0;
c_hi = A*k(1)^a;
i = 0;
dist_k = tol + 1;
dist_c = tol + 1;
while (i < I) & (dist_k > tol | dist_c > tol);
c(1) = (c_lo + c_hi)/2;
    for t = 1:T-1;
        k(t+1) = del_t*(A*k(t)^a - d*k(t) - c(t)) + k(t);
        c(t+1) = del_t*c(t)/s*(a*A*k(t)^(a-1) - p - d) + c(t);
        if c(t+1) > A*k(t+1)^a
           c(t+1) = A*k(t+1)^a;
        elseif c(t+1) < 0
               c(t+1) = 0;
        end
    end
    if k(T) > k_ss & c(T) < c_ss
        c_lo = c(1);
        elseif k(T) < k_ss & c(T) < c_ss
        c_hi = c(1);
        elseif k(T) < k_ss & c(T) > c_ss
        c_hi = c(1);
    end
dist_k = abs(k(T) - k(T-1));
dist_c = abs(c(T) - c(T-1));
i = i + 1;
end

% differential equations with y(1) = k, y(2) = c
dy = @(t,y) [A*y(1)^a - d*y(1) - y(2); ...     % dk/dt;
             y(2)/s*(a*y(1)^(a-1) - p - d)];   % dc/dt
% boundary conditions
bc = @(y0,yT) [y0(1) - k0;     % initial condition
               yT(2) - c_ss];  % terminal condition

solinit = bvpinit(linspace(0,T,10), [k0 c0]);   % initital guess
sol = bvp4c(dy, bc, solinit);                   % call solver
t = linspace(0,T)';  % time axis
y = deval(sol,t)';   % call solution values
k = y(:,1);    % transform variables
c = y(:,2);    % transform variables