% One_D_Parabolic_PDE.m
% Solution to Equation (6.18) for heat conduction (or diffusion). 
% This program calculates and plots the temperature (or concentration) 
% profiles, in a wall (or membrane) by solving the unsteady-state 
% equation using the function PARABOLIC1D.M.

clear;
redo = 1;
while redo
   clc;
   disp(' Solution of parabolic partial differential equation.')
   disp(' ')
   disp(' What type of problem are you solving?')
   disp('   1 - Heat conduction')
   disp('   2 - Diffusion')
   type = input(' Enter your choice : '); disp(' ')
   switch type
       case 1
       h = input(' Thickness of wall (m) = ');
       alpha = input(' Thermal diffusivity (m^2/s) = ');
       u0 = input(' Initial temperature of wall = ');
       y_label='Temperature'; TITLE =' Heat Conduction';
       case 2
       h = input(' Thickness of membrane (m) = ');
       alpha = input(' Diffusivity (m^2/s) = ');
       u0 = input(' Initial concentration in membrane = ');
       y_label='Concentration'; TITLE =' Diffusion';
   end
   disp(' ')
   tmax = input(' Maximum time of integration (s) = ');
   p = input(' Number of divisions in x-direction = ');
   q = input(' Number of divisions in t-direction = '); 
   disp(' '); disp(' Boundary conditions:')
   disp(' Condition at left side of wall (membrane):')
   disp('   1 - Dirichlet')
   disp('   2 - Neumann')
   disp('   3 - Robbins')
   bc(1,1) = input(' Enter your choice : ');
   if bc(1,1) < 3
      bc(1,2) = input(' Value = ');
   else
      disp(' u'' = (beta) + (gamma)*u')
      bc(1,2) = input(' Constant    (beta)  = ');
      bc(1,3) = input(' Coefficient (gamma) = ');
   end
   disp(' Condition at right side of wall (membrane):')
   disp('   1 - Dirichlet')
   disp('   2 - Neumann')
   disp('   3 - Robbins')
   bc(2,1) = input(' Enter your choice : ');
   if bc(2,1) < 3
      bc(2,2) = input(' Value = ');
   else
      disp(' u'' = (beta) + (gamma)*u')
      bc(2,2) = input(' Constant    (beta)  = ');
      bc(2,3) = input(' Coefficient (gamma) = ');
   end
   % Creating a vector of the initial condition
   U0 = [u0*ones(p+1,1)];
   % Calculating concentration profile
   [x,t,U] = parabolic1D(p,q,h/p,tmax/q,alpha,U0,bc);
   
   % Plotting K profiles
   K=10; count=fix(length(t)-1)/K;
   for i = 1:K
       new_U(:,i+1)=U(:,1+i*count);
   end
   clf; figure(1); plot(x,new_U);
   xlabel('Distance, m'); ylabel(y_label); title(TITLE);
   disp(' ')
   redo = input(' Repeat calculations (0/1) ? ');
end