function u=mgrid0(f)
[m,n]=size(f);N=m+1;
u=zeros(N+1,1);
if N==2
    u(2)=-f(1)/8;
    N
else
    h=1/N; u=0*ones(N+1,1);  om=0.66;
    r=( u(3:N+1)-2*u(2:N)+u(1:N-1) )/h^2 - f ; nr=norm(r);
    while nr > N*0.0002
     for i=1:5
         u(2:N) = u(2:N) +om*h^2*r/2;
         r=( u(3:N+1)-2*u(2:N)+u(1:N-1) )/h^2 - f ; nr=norm(r); nrv(i)=nr;
     end
     [N nrv],
     fgros=-r(3:2:N-1);
     ucorr=mgrid0(fgros);
     u(3:2:N-1)=u(3:2:N-1)+ucorr(2:N/2);
     u(2:2:N)=u(2:2:N)+0.5*(ucorr(1:N/2)+ucorr(2:N/2+1) );
     r=( u(3:N+1)-2*u(2:N)+u(1:N-1) )/h^2 - f ; nr=norm(r);
    end
end
