function err=usediff2(J);
% USEDIFF2  Function which demonstrates use of DIFFUSION2 on manufactured
% example where exact solution is known:
%     div( a(x,y) grad u) + f = 0
% where  u(x,y)=exp(y)-x   and  a(x,y)=1+0.75 sin(3*x).   Returns maximum 
% error on grid of size J.
% Example:
%   >> usediff2(40), usediff2(80), usediff2(120), usediff2(160)
% ELB 4/24/05; update which uses DIFFUSION2 and a nonsymmetical
% manufactured example.
% See also: DIFFUSION2

dx=1/J;  dy=dx;  xx=0:dx:1;  yy=xx;  % xx,yy have length J+1
[XX,YY]=ndgrid(xx,yy);
fmatrix=f(XX(2:J,2:J),YY(2:J,2:J));  % eval f at interior points
[XL,YL]=ndgrid(dx/2:dx:1-(dx/2),yy(2:J));  % left staggered grid
[XD,YD]=ndgrid(xx(2:J),dy/2:dy:1-(dy/2));  % down staggered grid
aLEFT=a(XL,YL);  aDOWN=a(XD,YD);
U=diffusion2(J,aLEFT,aDOWN,fmatrix,@g);  % pass function *handle* for g
uex=g(XX,YY);  err=max(max(abs(U-uex)));  % exact solution and error

function z=a(x,y);  z=1+0.75*sin(3*x);  % diffusivity

function z=f(x,y);  z= 2.25*cos(3*x) - (1+0.75*sin(3*x)) .* exp(y);  % source

function u=g(x,y);  u=exp(y)-x;   % both boundary value and exact solution