% Initialize b, x, jac1, jac2, iter to demonstrate testfmg
% Written by J. Demmel, UC Berkeley, July 10, 1993
figure(1)
hold off, clg
k=4;
kk=2^k+1;
[x,y]=meshdom(1:kk,1:kk);
b=sin((x-1)*3*pi/kk).*sin((y-1)*pi/kk);
b=rand(size(b))-.5;
b=b*0;
p1=((kk-1)/8:3*(kk-1)/8)+1;
p2=(5*(kk-1)/8:7*(kk-1)/8)-1;
b(p1,p2)=b(p1,p2) -2*ones(prod(size(p1)),prod(size(p2)));
b(p2,p1)=b(p2,p1) + ones(prod(size(p2)),prod(size(p1)));
[nb,mb]=size(b);
b(:,1)=zeros(nb,1);b(:,mb)=zeros(nb,1);
b(1,:)=zeros(1,mb);b(nb,:)=zeros(1,mb);
xtrue = b;
b=zeros(size(b));
b(2:kk-1,2:kk-1) = 4*xtrue(2:kk-1,2:kk-1) ...
-xtrue(1:kk-2,2:kk-1) ...
-xtrue(3:kk ,2:kk-1) ...
-xtrue(2:kk-1,1:kk-2) ...
-xtrue(2:kk-1,3:kk );
hold off; clg
% x=b;b=zeros(size(b));
subplot(1,2,1); mesh(xtrue), title('true solution')
subplot(1,2,2); mesh(b), title('right hand side')
x=zeros(size(b));
jac1=1;jac2=1;iter=20;