function resid = mydpresid(model,fspace,s,v,x,c,scoord_resid)

% SET PARAMETER & SHOCK DISTRIBUTION DEFAULTS
  tol       = optget('dpsolve','tol',sqrt(eps));
  maxit     = optget('dpsolve','maxit',500);
  nres      = optget('dpsolve','nres',10);
  showiters = optget('dpsolve','showiters',1);
  algorithm = optget('dpsolve','algorithm','newton');

  if ~isfield(model,'e'); model.e=0; e=0;
    else, e=model.e; end;
  if ~isfield(model,'w'); model.w=1; delw=model.discount;
    else, delw=model.discount*model.w; end;
  if isfield(model,'horizon') & model.horizon<inf, algorithm='finite'; end
  if isfield(model,'actions'), x=model.actions; end
  if isfield(model,'explicit'), explicit=model.explicit; else, explicit=0; end
  func=model.func;
  params=model.params;
  if nargin<3, s=gridmake(funnode(fspace)); end
  if nargin<4 | isempty(v), v = zeros(size(s,1),1); end

% DETERMINE NUMBER OF DIMENSIONS & COORDINATES
  n  = fspace.n;                % number of collocation coordinates by state dimension 
  ns = prod(n);                 % number of collocation states
  ds = length(n);               % dimension of state space 
  dx = size(x,ndims(x));               % dimension of action space
  
% COMPUTE COLLOCATION NODES AND INTERPOLATION MATRIX 
  PhiS = funbasx(fspace,s);          % collocation matrix
  [F,G,X] = GetFG(s,fspace,model);   % basis functions for discrete states
   

% COMPUTE RESIDUAL
s = gridmake(scoord_resid);
if isempty(F)
  if explicit
    x=[];
  else
    x = funeval(funfitxy(fspace,PhiS,x),fspace,scoord_resid);   % rough guess for actions at evaluation points
  end
  [v,x] = vmax(s,x,c,fspace,model,F,G);                 % values and actions at evaluation points
else
  [ss,xx]=gridmake(s,X);
  v = valfunc(c,fspace,ss,xx,e,delw,func,params);
  [v,j] = max(reshape(v,size(s,1),size(X,1)),[],2); % values and actions at evaluation points
  x = X(j,:);
end
resid = v-funeval(c,fspace,scoord_resid);   % residual at evaluation points
n_resid = [];
for ii = 1:length(scoord_resid)
    n_resid(ii) = length(scoord_resid{ii});
end
resid = reshape(resid,[n_resid 1]);         % reshape residual for plotting

  

    
% VMAX - Solves Bellman Equation
  function [v,x,vc] = vmax(s,x,c,fspace,model,F,G,X)
  % CONTINUOUS ACTIONS
  if isempty(F)
    % SET CONVERGENCE PARAMETER DEFAULTS
    tol           = optget('dpsolve_vmax','tol',5e-8);
    maxit         = optget('dpsolve_vmax','maxit',50);
    maxbacksteps  = optget('dpsolve_vmax','maxbacksteps',0);
    lcpmethod     = optget('dpsolve_vmax','lcpmethod','minmax');

    e=model.e;
    delw = model.discount*model.w;
    func=model.func;
    params=model.params;
    
    if isempty(x) % explict solution to control available
      dv=squeeze(funeval(c,fspace,s,eye(size(s,2))));
      x=feval(func,'x',s,dv,[],params{:});
      v = valfunc(c,fspace,s,x,e,delw,func,params);
    else
      % COMPUTE BOUNDS      
      [xl,xu] = feval(func,'b',s,x,[],params{:}); 

      % SOLVE FIRST ORDER CONDITIONS
      for it=1:maxit
        [v,vx,vxx] = valfunc(c,fspace,s,x,e,delw,func,params);
        [vx,deltax] = lcpstep(lcpmethod,x,xl,xu,vx,vxx);
        err = max(abs(vx),[],2);
        if all(err<tol), break, end;
        eold = inf;
        if maxbacksteps<1
          x = x+deltax;
        else
          for k=1:maxbacksteps
            xnew = x + deltax;
            [v,vx] = valfunc(c,fspace,s,xnew,e,delw,func,params);
            vx = lcpstep(lcpmethod,x,xl,xu,vx);
            enew = max(abs(vx),[],2);
            ind = find(eold>enew & enew>err);
            if isempty(ind), break; end
            eold = enew;
            deltax(ind,:) = deltax(ind,:)/2;
          end 
          x = xnew;
        end
      end
    end

    % COMPUTE dv/dc 
    if nargout>2
      ns = size(s,1);
      g = feval(func,'g',s,x,e(ones(ns,1),:),params{:});
      vc = delw(1)*funbas(fspace,g);
      for k=2:length(delw)
        g = feval(func,'g',s,x,e(k+zeros(ns,1),:),params{:});
        vc = vc + delw(k)*funbas(fspace,g);
      end
    end
  % DISCRETE ACTIONS
  else
    [ns,nc] = size(G);
    nx = length(model.actions);
    ns = ns/nx;
    [v,j] = max(reshape(F+G*c,ns,nx),[],2);
    x = X(j,:);
    if nargout>2
      i = (j-1)*ns + (1:ns)';   
      vc = reshape(G(i,:),ns,nc);
    end
  end
  
 
% VALFUNC  Evaluates Bellman Optimand
  function [v,vx,vxx]=valfunc(c,fspace,s,x,e,delw,func,params);
  
  % COMPUTE LOCAL CONSTANTS
    [ns,ds] = size(s);
    dx  = size(x,2);
    dxx = dx*dx;
    K = length(delw);

    if nargout<2
      v = feval(func,'f',s,x,[],params{:});
      for k=1:K
        kk = k + zeros(ns,1); 
        g = feval(func,'g',s,x,e(kk,:),params{:});
        v  = v  + delw(k)*funeval(c,fspace,g);
      end
    elseif nargout<3
      [v,vx] = feval(func,'f',s,x,[],params{:});
      vx = reshape(vx,ns,1,dx);
      for k=1:K
        kk = k + zeros(ns,1); 
        [g,gx] = feval(func,'g',s,x,e(kk,:),params{:});
        [vnext,vs] = fund(c,fspace,g,1);
        v  = v  + delw(k)*vnext;
        vx = vx + delw(k)*arraymult(vs,gx,ns,1,ds,dx);
      end
      clear g gx 
      vx = reshape(vx,ns,dx);
    else
      [v,vx,vxx] = feval(func,'f',s,x,[],params{:});
      vx  = reshape(vx,ns,1,dx);
      vxx = reshape(vxx,ns,dx,dx);
      for k=1:K
        kk = k + zeros(ns,1);
        [g,gx,gxx] = feval(func,'g',s,x,e(kk,:),params{:});
        [vnext,vs,vss] = fund(c,fspace,g,1);
        v   = v   + delw(k)*vnext;
        vx  = vx  + delw(k)*arraymult(vs,gx,ns,1,ds,dx);
        vxx = vxx + delw(k)*(reshape(arraymult(vs,gxx,ns,1,ds,dxx),ns,dx,dx) ...
          + arraymult(permute(gx,[1 3 2]),arraymult(vss,gx,ns,ds,ds,dx),ns,dx,ds,dx));
      end  
      clear g gx gxx vss
      vx  = reshape(vx,ns,dx);
      vxx = reshape(vxx,ns,dxx);
    end

  
% GETFG - Computes Reward and Discounted Expected Basis for Discrete Choice Models
  function [F,G,X] = GetFG(s,fspace,model)
  
  if ~isfield(model,'actions'), F=[];G=[]; X=[]; return; end

  e = model.e;
  delw = model.discount*model.w;
  func=model.func;
  params=model.params;

  K = length(delw);
  X = gridmake(model.actions);
  [ss,xx] = gridmake(s,X);
  ns=size(ss,1);
  
  F = feval(func,'f',ss,xx,[],params{:});
  g = feval(func,'g',ss,xx,e(ones(ns,1),:),params{:});
  G=delw(1)*funbas(fspace,g);
  for k=2:K
    g = feval(func,'g',ss,xx,e(k + zeros(ns,1),:),params{:});
    G = G + delw(k)*funbas(fspace,g);
  end    