%Homework #2
%Creator		Kyaw Tun
%varsion		1			original code
function main1()

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Set system parameter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global a b NSTEP TOLERANCE MAX_ITERATION E I delta_x;
a = 0;
b = 4;
NSTEP = 41;
TOLERANCE = 10^(-10);
MAX_ITERATION = 20;
global u u_old;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Initialization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
I = 1;
E = 0; E_old = 0;
delta_x = (b - a) / NSTEP;
for i = 1: NSTEP;
   u_old(i) = 1;
   x(i) = a + (i-1) * delta_x;
   u(i) = 2;
end


while I < MAX_ITERATION			%Main loop start
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Calculate u(i) 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i = 1: NSTEP
   u(i) = exp(-Integrator(i));
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Check tolerance
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
E = 0;
for i = 1: NSTEP
   E = E + abs(u(i) - u_old(i));	
end
if E <= TOLERANCE
   Reporter(1);
   return;
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Registertor Office
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
I = I + 1;
if E_old > E
   fprintf('>');
elseif   E_old < E
   fprintf('<');
else
   fprintf('-');
end
E_old = E;
for i = 1: NSTEP
   if isnan(u(i))
      Reporter(-1, 'Not a number');
      return;
   elseif isinf(u(i))
      Reporter(-1, 'u(i) is infinate');
      return;
   end
   u_old(i) = u(i);
end
if I >= MAX_ITERATION - 1
   fprintf('\nMaxinum number of iteration %d is reach.', MAX_ITERATION);
   y = input('\nEnter next number of iterations (<100) or press return to exit ');
   if isnumeric(y) & (y > 0) & (y < 100) 
      MAX_ITERATION = MAX_ITERATION + y;
   else
      return;
   end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%% Registertor Office end



end										%main loop end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%  Utility functions %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Reporter(x,s)
%this function terminate or reset the program 
global E I NSTEP u a delta_x;
if x == 1
   fprintf('\nResult obtaine!\n within %d iteration.\n Error is %e.\n', I, E);
    [filename, pathname] = uiputfile('*.doc', 'Result file name');
    FileName = [pathname, filename, '.doc'];
   [fid msg] = fopen(FileName, 'w');
   if fid == -1
      disp(msg)
      disp(['Result cannot be writed in ', FileName, 'It may be opened.']);
      return
   else  
      fprintf('\nResult is writed down in %s\n', FileName);
      fprintf(fid, '\n Results\n Error = %e\n Number of iteration = %d\n', E, I);
	   for i = 1: NSTEP
         fprintf(fid, '\n u(%d) = %e', i, u(i));
         X(i) = (i-1) * delta_x + a;
      end 
      fclose(fid);  
   end 
   plot(X, u);
   xlabel('x');
   ylabel('u(x)');
   title('plot of u');
   return;
end
if x == -1
   disp(s);
end
if x == 2
   disp(s);
   return;
end



function y = Integrator(x)
%Exsecute integration using trapezoidal rule
global a b NSTEP delta_x;
Integrand_sum = 0;
x1 = a + delta_x * x;
for i = 2: (NSTEP - 1)
   x2 = (i - 1) * delta_x;
   Integrand_sum = Integrand_sum + Integrand(x1, x2, i);   
end
y = delta_x / 2 * (Integrand(x1, a, 1) + 2 * Integrand_sum + Integrand(x1, b, NSTEP));



function y = Integrand(x1, x2, i)
global u_old;
%Evalute integrand
y = 1/2 * exp( - abs(x1-x2)) * u_old(i);