function err = errfun(x,D2,pfo,Re,D,sTR2,cotT,R)
% ERRFUN - Navier-Stokes error function.
%
% Syntax: err = errfun(x,D2,pfo)
%
% This function returns the error in the Navier-Stokes equations, so 
% to solve the N-S equations means to find x such that errfun(x,...) = 0.

% Recover grid size
[M,N] = size(D2);

% Turn x from a vector into a matrix
x = reshape(x,2*M,N);

% Pick out the x1=psi and the x2=D2(psi) components
x1 = x(1:M,:);
x2 = x(M+(1:M),:);

% The equation x2 = D2(x1)
f1 = D2 * x1 - x2;

% Boundary conditions: psi = 0 at inner boudnary and symmetry lines, free flow at outer.
f1(1,:) = x1(1,:);
f1(:,1) = x1(:,1);
f1(:,end) = x1(:,end);
f1(end,:) = x1(end,:) - pfo;

% The equation D2(D2(psi)) = 0 (only works for Re = 0 at this point)
A = D2 - Re./sTR2*((D.dth*x1)*D.dr - (D.dr*x1) * D.dth + ( 2* cotT.* (D.dr*x1) - 2*(D.dth*x1)./R) * eye(D.dr));
f2 = A * x2;

% Boundary conditions: D2(psi) = 0 at outer boundary and symmetry lines; at inner
% boundary we have D.dr(psi) = 0
f2(1,:) = x1(1,:) - x1(2,:);
f2(:,1) = x2(:,1);
f2(:,end) = x2(:,end);
f2(end,:) = x2(end,:);

% Stack up the two error components and return as a vector.
err = [f1;f2];
err=err(:);

return