/////////////////////
// TP1, exercice 4 //
/////////////////////


//Le nombre de pas de discrétisation
N=30;

//Le pas d'espace
dx=1/(N+1);


////////////////////////////////////////////////////


//Définition de la matrice du problème

//Definition des blocs
bloc_1=-4*diag(ones(1,N)) + diag(ones(1,N-1),1) + diag(ones(1,N-1),-1);
bloc_2=diag(ones(1,N));

//On définit une matrice vide que l'on va remplir
D=zeros(N^2,N^2);

//On commence par remplir les blocs diagonaux
for i=0:N-1
  D(1+i*N:(i+1)*N , 1+i*N:(i+1)*N)=bloc_1;
end;
//Puis on remplit les blocs sur- et sous-diagonaux
for i=1:N-1
  D([1+i*N : (i+1)*N]  ,  [1+(i-1)*N : i*N])=bloc_2;
  D([1+(i-1)*N : i*N]  ,  [1+i*N : (i+1)*N])=bloc_2;
end;
D=D/(dx^2);

//La matrice de -(Laplacien)+20*(Identité)
A=-D+20*diag(ones(1,N^2));

//Définition du second membre
f=zeros(N^2,1);
for i=1:N
  f(i)=f(i)+i*dx;
  f(i*N)=f(i*N)+1;
  f(N^2-N+i)=f(N^2-N+i)+i*dx;
end;
f=f/(dx^2);


///////////////////////////////////////////////


//On resoud le schéma numérique
sol_vect=A\f;

//On écrit la solution sous forme de matrice (sol_vect était un vecteur colonne)
sol_mat=zeros(N+2,N+2);
sol_mat(1,N+2)=1;
sol_mat(N+2,N+2)=1;
for i=1:N
  sol_mat(1,i+1)=i*dx;
  sol_mat(N+2,i+1)=i*dx;
  sol_mat(i+1,N+2)=1;
  for j=1:N
	sol_mat(i+1,j+1)=sol_vect(j+N*(i-1));
  end;
end;


//////////////////////////////////////////////


//Tracé de la solution approchée
fenetre=scf(0);clf();//Une nouvelle fenêtre
fenetre.color_map=autumncolormap(50);//Nouveau choix de couleurs pour l'affichage
grayplot(dx*[0:N+1],dx*[0:N+1],sol_mat);

fenetre=scf(1);clf();//Une nouvelle fenêtre
fenetre.color_map=autumncolormap(50);//Nouveau choix de couleurs pour l'affichage
plot3d1(dx*[0:N+1],dx*[0:N+1],sol_mat);
