%
% MATLAB for DUMMIES 11-12
% Exercice 8
%
% Solution detaillee
%  Vincent Legat, Damien Bouvy, Denis-Gabriel Caprace, Gilles De Neyer, 
%  Kevin Degraux, Mohammad Iranmanesh, Christophe Lorant, Remy Moulin, 
%  Matthew Philippe, Karim Slaoui, Olivier Thiry, Laurent Van Eesbeeck,
%  Thomas Vanderstraeten, Antoine Weynants
%
% -------------------------------------------------------------------------
%
%    Ceci est une solution SIMPLE mais parfaitement correcte :-)
%    Elle n'est pas la plus rapide, mais pas si lente que cela...
%
%    Le resultat du concours sera disponible pour l'An Neuf...
%    Certains etudiants ont ete tres ingenieux.
%    Les meilleurs programmes seront mis sur le site du cours.
%
function [Us,Umax,Xmax] = theBigFish(alpha);

n = 4*alpha+1;
h = 2/(n-1);
big = n*n*n;
%
% -1- Allocation d'une matrice CREUSE de taille n^3
%     C'est l'ingredient essentiel pour l'efficacite !
%
A = sparse(big,big);
b = zeros(big,1);
for i=1:big
    A(i,i) = 1;
end
%
% -2- Ecriture des equations (on pourrait optimiser ici :-)
%     Mais, c'est tellement plus simple ainsi
%
for i=2:n-1
    for j =2:n-1
        for k =2:n-1
            index = n*n*(k-1)+n*(j-1)+i;
            A(index,index)     = -6; 
            A(index,index+1)   = 1;
            A(index,index-1)   = 1;
            A(index,index+n)   = 1;
            A(index,index-n)   = 1;
            A(index,index+n*n) = 1;
            A(index,index-n*n) = 1;
            if (i > n/2 && j > n/2 && k > n/2) b(index) = -10;  end
        end      
    end
end
A = A/(h^2); 
%
% -3- Resolution du systeme creux
%     Pas evident de faire mieux que le solveur creux de MATLAB
%     Quoiqu'un solveur multigrille ou des gradients conjuges peuvent
%     etre efficaces ici (malgre la taille reduite du probleme)
%
U = A\b;
Us = reshape(U,n,n,n);
%
% -4- Calcul de la position du maximum
%     Qui ne se trouve pas en 0.5,0.5,0.5 !
%     Le tout petit piege du probleme...
%
[Umax, Index] = max(U);
X = h*(mod(Index,n)-1) - 1; 
Xmax = [X X X];

end