function multigrid()
% Solving Equation grad * (a(x,y) * grad(u)) = f(x,y)
p = 10;
n = 2^p + 1; % the finest mesh is n-by-n
n1 = n - 1;
nsmooth = 2; % # of iterations for smoothing
mytol = 1.0e-9;
hmax = 1/2; % the maximal mesh step
hmin = 1/n1; % the minimal mesh step

t = linspace(0,1,n);
[x y] = meshgrid(t,t);
% set f(x,y) so that the exact solution is the given one 
uexact = x.*(1-x).*y.*(1-y); % u exact
a = 1 + x + 2*y.^2;
f = zeros(n);
f(2:n1,2:n1) = diff_oper(hmin,uexact,a);
r = f; % the residual
%% multigrid - Gauss Seidel V- cycle

u = zeros(n);  % the initial guess
res(1) = max(max(abs(diff_oper(hmin,u,a) - f(2:n1,2:n1))));
err(1) = max(max(abs(u - uexact)));
iter = 0;
fprintf('V-cycle #%d: residual = %d, error = %d \n',iter,res(iter + 1),err(iter + 1));

tic;
while res(iter + 1) > mytol && iter < 40
    e = zeros(n); % the initial guess for the error
    % smoothing and restriction
    h = hmin;
    step = 1;
    while h < hmax - eps  % restriction
        h = h * 2;
        step = step * 2;
    end

    % solution on the 3-by-3 grid
    m = 1 + step; % the index of the middle point
    h = hmax;
    e(m,m) = r(m,m)*2*h^2/(a(1,m) + a(m,1) + 4*a(m,m) + a(m,n) + a(n,m));

    % prolongation and post-smoothing
    I = [1 : step : n];
    while h > hmin + eps
        h = h/2;
        step = step/2;
        % prolongation by linear interpolation
        Ifine = [1 : step : n];
        e(Ifine,Ifine) = interp2(x(I,I),y(I,I),e(I,I),x(Ifine,Ifine),y(Ifine,Ifine));
        % post-smoothing
        I = Ifine;
        e(I,I) = gauss_seidel(nsmooth,h,e(I,I),a(I,I),r(I,I));
    end
    iter = iter+1;
    u = u + e; % correct the solution
    r(2:n1,2:n1) = f(2:n1,2:n1) - diff_oper(hmin,u,a); %r(2:n1,2:n1) - diff_oper(h,e,a); % the residual
    res(iter + 1) = max(max(abs(r)));
    err(iter + 1) = max(max(abs(u - uexact)));
    fprintf('iter #%d: residual = %d, error = %d \n',iter,res(iter + 1),err(iter + 1));
end
tottime = toc;
fprintf('Gauss-Seidel: V-cycles = %d, residual = %d, CPUtime = %d\n',iter,res(iter + 1),tottime);
figure;
clf;
hold on;
plot([0:iter],res,'Linewidth',1,'marker','*','color','b');
set(gca,'YScale','log','Fontsize',24)
title('Gauss-Seidel','FontSize',24)
xlabel('V cycle #','Fontsize',20);
ylabel('Residual','FontSize',20);
grid
hold off


end
%%
function u=gauss_seidel(niter,h,u,a,f);
[m n]=size(u);
% coefficients for the finite difference scheme
% ((au_x)_x + (au_y)_y + f)(i,j) = 
% (1/(2h^2))(amx*u(i-1,j) + apx*u(i+1,i) + amy*u(i,j-1)+apy*u(i,j+1) -
%            acxy*u(i,j)) + f(i,j) = 0

amx=zeros(m,n);
apx=zeros(m,n);
amy=zeros(m,n);
apy=zeros(m,n);
acxy=zeros(m,n);
amx(2:m,1:n)=a(1:m-1,1:n)+a(2:m,1:n);
apx(1:m-1,1:n)=amx(2:m,1:n);
amy(1:m,2:n)=a(1:m,1:n-1)+a(1:n,2:n);
apy(1:m,1:n-1)=amy(1:m,2:n);
acxy(2:m-1,2:n-1)=amx(2:m-1,2:n-1)+apx(2:m-1,2:n-1)+amy(2:m-1,2:n-1)+apy(2:m-1,2:n-1);

for k=1:niter
    for i=2:m-1
        for j=2:n-1
            u(i,j)=(f(i,j)*2*h^2+u(i-1,j)*amx(i,j)+u(i+1,j)*apx(i,j)+...
                u(i,j-1)*amy(i,j)+u(i,j+1)*apy(i,j))/acxy(i,j);
        end
    end
end
end


%%
function Lu=diff_oper(h,u,a)
[m n]=size(u);

amx=zeros(m,n);
apx=zeros(m,n);
amy=zeros(m,n);
apy=zeros(m,n);
acxy=zeros(m,n);
amx(2:m,1:n)=a(1:m-1,1:n)+a(2:m,1:n);
apx(1:m-1,1:n)=amx(2:m,1:n);
amy(1:m,2:n)=a(1:m,1:n-1)+a(1:n,2:n);
apy(1:m,1:n-1)=amy(1:m,2:n);
acxy(2:m-1,2:n-1)=amx(2:m-1,2:n-1)+apx(2:m-1,2:n-1)+amy(2:m-1,2:n-1)+apy(2:m-1,2:n-1);

Lu=(acxy(2:m-1,2:n-1).*u(2:m-1,2:n-1)-amx(2:m-1,2:n-1).*u(1:m-2,2:n-1)...
    -apx(2:m-1,2:n-1).*u(3:m,2:n-1)-amy(2:m-1,2:n-1).*u(2:m-1,1:n-2)...
    -apy(2:m-1,2:n-1).*u(2:m-1,3:n))/(2*h^2);

end


