format short e; format compact;
Lx=1;Ly=1;
dx=Lx/(n);
dy=dx;
tol=dx^3;
x=0:dx:1;
y=0:dy:1;
[X,Y]=meshgrid(x,y);
f=-sin(pi*X).*sin(pi*Y);
vex=-f/pi^2/2;
v=0*vex;
uancien=v;
% iterations de GS red/black
errgsrb=100;
itegsrb=0;
t0=cputime;
while(errgsrb>tol)
v0=v;
itegsrb=itegsrb+1;
indix=2:2:n;
indiy=2:2:n;
t=v(indix+1,indiy)+v(indix-1,indiy)+v(indix,indiy+1)+v(indix,indiy-1);
v(indix,indiy)=(t-dx^2*f(indix,indiy))/4;
indix=3:2:n-1;
indiy=3:2:n-1;
t=v(indix+1,indiy)+v(indix-1,indiy)+v(indix,indiy+1)+v(indix,indiy-1);
v(indix,indiy)=(t-dx^2*f(indix,indiy))/4;
indix=2:2:n;
indiy=3:2:n-1;
t=v(indix+1,indiy)+v(indix-1,indiy)+v(indix,indiy+1)+v(indix,indiy-1);
v(indix,indiy)=(t-dx^2*f(indix,indiy))/4;
indix=3:2:n-1;
indiy=2:2:n;
t=v(indix+1,indiy)+v(indix-1,indiy)+v(indix,indiy+1)+v(indix,indiy-1);
v(indix,indiy)=(t-dx^2*f(indix,indiy))/4;
errgsrb=max(max(abs(v-v0)))/max(max(abs(v)));
ax=[0,1,0,1,0,1/pi^2];
%surf(x,y,v);axis(ax);pause(0.01);
end;
t1=cputime-t0;
disp([n,itegsrb,t1]);




