% Example6_1.m
% Solution to Example 6-1. This program calculates and plots
% the solution of the Laplace (or Poisson) equation within a 
% rectangular plate by finite difference method, using the 
% function ELLIPTIC.M.

clf
redo = 1;
while redo
   clear; clc
   bcdialog = [' Lower x boundary condition:'
               ' Upper x boundary condition:'
               ' Lower y boundary condition:'
               ' Upper y boundary condition:'];           
   disp('         Solution of the Laplace (or Poisson) equation')
   disp('    (Two-dimensional elliptic partial differential equation)')   
   disp(' '); disp(' ')
   disp('                          Upper y boundary          '          )
   disp('              ______________________________________'          )
   disp('             |                                      |'         )   
   disp('             |                                      |'         )   
   disp('             |                                      |'         )
   disp('             |                                      |'         )
   disp('    Lower x  |                                      | Upper x' )
   disp('    boundary |                                      | boundary')
   disp('             |                                      |'         ) 
   disp('             |                                      |'         )
   disp('             y                                      |'         )   
   disp('             |___x____________length________________|'         )
   disp('            0                                        '         )
   disp('                          Lower y boundary           '         )

   disp(' ')
   disp(' ')
   length = input(' Length of the plate (x-direction) (m) = ');
   width = input(' Width of the plate  (y-direction) (m) = ');
   p = input(' Number of divisions in x-direction = ');
   q = input(' Number of divisions in y-direction = ');
   f = input(' Right hand side of the equation (f) = ');
   disp(' '), disp(' ')
   disp(' Boundary conditions:')
   for k = 1:4
      disp(' ')
      disp(bcdialog(k,:))
      disp(' 1 - Dirichlet')
      disp(' 2 - Neumann')
      disp(' 3 - Robbins')
      bc(k,1) = input(' Enter your choice : ');
      if bc(k,1) < 3
         bc(k,2) = input(' Value = ');
      else
         disp(' u'' = (beta) + (gamma)*u')
         bc(k,2) = input(' Constant    (beta)  = ');
         bc(k,3) = input(' Coefficient (gamma) = ');
      end
   end
   [x,y,T] = elliptic(p,q,length/p,width/q,bc,f);  
   figure(1);
   surf(x,y,T)
   xlabel('x (m)')
   ylabel('y (m)')
   zlabel('T (deg C)')
   colorbar
   shading interp
   figure(2);
   pcolor(x,y,T)
   xlabel('x (m)')
   ylabel('y (m)')
   zlabel('T (deg C)')
   view(2)   
   shading interp;
   colorbar
   axis('equal')
   axis('tight')

   print_matrix = input(' Do you want to print the matrix of results (0/1) ? ');
   % Flip the matrix to arrange the position of values
   % to correspond to the points on the grid
   for i=1:q+1
   TN(i,:)=T(q+2-i,:);
   end 
   if print_matrix; TN, end
   disp(' ')
   redo = input(' Repeat calculations (0/1) ? ');
end