function [interp,exitflag,pstar] = SolveStorageDL(model,interp,options)
% SOLVESTORAGEDL Solves the storage model by fixed point as proposed by Deaton and Laroque

%% Initialization
params              = num2cell(model.params);
[a, b, delta, k, r] = params{:};
beta                = (1-delta)/(1+r);

e = model.shocks.e;
w = model.shocks.w;

demand       = @(p) (p-a)/b;
invdemand    = @(d) a+b*d;

A            = interp.s;

PA           = invdemand(A); % Precalculation for speed

Display      = strcmpi(options.Display,'iter');
InterpMethod = options.InterpMethod;
MaxIter      = options.MaxIter;
TolX         = options.TolX;

n            = length(A);
kshocks      = size(model.shocks.e,1);

% First guess
if isfield(interp,'cx')
  PriceInterp = interp.cx{2};
else
  PriceInterp = griddedInterpolant(A,max(PA,0),InterpMethod);
end

if isfield(interp,'x')
  P = interp.x(:,2);
else
  P = max(PriceInterp(A),PA);
end

dis      = inf;
exitflag = 1;
Iter     = 0;

%% Successive iterations
if Display
  fprintf(1,'Successive approximation\n');
  fprintf(1,' Iteration  Residual\n');
end

while(dis > TolX && Iter < MaxIter)
  Iter = Iter + 1;
  Pold = P;

  %% Calculate next-period availability
  S     = max(A-demand(P),0);
  Anext = reshape(bsxfun(@plus,(1-delta)*S',e),kshocks*n,1);

  %% Update the price and its approximation
  Pnext = max(PriceInterp(Anext),invdemand(Anext));
  P     = max(PA,beta*reshape(w'*reshape(Pnext,kshocks,n),n,1)-k);
  PriceInterp.Values = P;

  dis = norm(P-Pold);
  if Display, fprintf(1,'%9i   %7.2E\n',Iter,dis); end
end

S                  = max(A-demand(P),0);
StockInterp        = griddedInterpolant(A,S,InterpMethod);
interp.x           = [S P];
interp.cx          = {StockInterp; PriceInterp};

if Iter==MaxIter
  warning('Storage:FailureDLSolve','Failure to solve for REE, residuals= %4.2e',dis)
  exitflag = 0;
end

if nargout==3
  pstar = w'*max(PriceInterp(e),invdemand(e))*(1-delta)/(1+r)-k;
end