%                     Lab Assignment #1
%  
%   Program the following Matlab code. The best way is to make a 
%   file (filename: lab1.m) of the commands as they are shown below.  
%   Startup matlab on your system and at the matlab prompt 
%   (typically >) type: lab1
%   The program should start up and prompt you for input.
%
%   The system you are solving is 
%                n+1    n          n       n      n
%               u    = u  - CFL(a u   + b u  + c u   )
%                j      j          j-1     j      j+1
%   with periodic boundary conditions
%
%   It is the simple wave equation on a periodic domain.  We have chosen a 
%   three point differencing in space and explicit Euler time differencing.
%             
%   Problem #1:  For nmax = 10, CFL = 1.0, run the program with
%      a)  a = -0.5, b = 0.0, c = 0.5 and describe what happens.
%          Short descriptions please.
%          Hand sketch the 5th and 10th curve plot of u vrs x
%          or if you have printer capability print the results.
%      b)  Repeat a) for a = -1, b = 1, c = 0
%      c)  Repeat a) for CFL = 0.7
%      d)  Repeat a) for CFL = 1.4
%
%   Problem #2:  For nmax = 10, a = -1, b = 1, c = 0  run the program with
%      a)  CFL = 0.7 and describe what happens.
%          Hand sketch the 5th and 10th curve plot of u vrs x
%          or if you have printer capability print the results.
%      b)  Repeat a) CFL = 1.4
%        
%   Problem #3: 
%          For nmax = 10, CFL=1.0, run the program with
%      a)  a = -0.9, b = 0.8, c = 0.1 and describe what happens.
%          Short descriptions please.
%          Hand sketch the 5th and 10th curve plot of u vrs x
%          or if you have printer capability print the results.
%      b)  Repeat a) for a = -0.99, b = 0.98, c = 0.01

help lab1;

%% Start of program ::::::::::::::::::::::::::::::::::::::::::::

jmax = 51;       

% The grid (simply the points 1 to jmax)
x = [1:jmax];

%  Initial Conditions
u1  = zeros(jmax,1);   
un1 = zeros(jmax,1);   
u2  = zeros(jmax,1);   
un2 = zeros(jmax,1);   

%  set spike at middle
u1((jmax-1)/2) = 1.0;
u2((jmax-1)/2) = 1.0;

%  Input
nmax = 10; 
%nmax = input('Enter the number of steps to take :: ');

CFL1 = 1.4;  
CFL2 = .7;
%CFL = input('Enter CFL :: ');

a1 = -0.5;  
b1 = 0.; 
c1 = 0.5;
a2 = -0.5;  
b2 = 0.; 
c2 = 0.5;

%fprintf('Enter a,b,c of B(a,b,c) \n');
%a = input('Enter of  a :: ');
%b = input('Enter of  b :: ');
%c = input('Enter of  c :: ');
clf;
%  iterate nmax times 
for iter = 1:nmax 
  
  % interior scheme
  %  The points we are updating (integrating) are j = 2 to j = jmax-1
  %  The boundary points j=1 and j=jmax are overloaded 
  %  with j=jmax-1 and j=2 respectively to simulate periodicty, i.e. to
  %  elliminate end effects

  j=[2:jmax-1]; % the interior index array

  %  load new values into array un
  un1(j) = u1(j) - CFL1 * (a1*u1(j-1) + b1*u1(j) + c1*u1(j+1));
  un2(j) = u2(j) - CFL2 * (a2*u2(j-1) + b2*u2(j) + c2*u2(j+1));
  % Note the above two lines are equivalent to
  %    for j=2:jmax-1
  %       un(j) = u(j) - CFL * (a*u(j-1) + b*u(j) + c*u(j+1));
  %    end
  % But, using the vector form is faster (much!)
  
  % boundary scheme enforces periodicity
  un1(1)    = un1(jmax-1);  
  un1(jmax) = un1(2);
  un2(1)    = un2(jmax-1);  
  un2(jmax) = un2(2);
  
  % update u
  u1        = un1;
  u2        = un2;
  %  If you have plotting capability use the next statement 
  subplot(2,1,1);
  plot(x,u1);
  title(strcat('Wave Equation at n=',num2str(iter)));
  xlabel('x');ylabel('u');
  axis([1 jmax min(u1) max(max(u1),1)]);
  legend('U1','U22');
  subplot(2,1,2);
  plot(x,u2);
  axis([1 jmax min(u2) max(max(u2),1)]);
  pause;
  %  Note!! the scales on the plot may be changing dynamically
  
  %  The next statement outputs data to the screen, the lack 
  %  of a trailing ";" is correct, causing the output.
  u1'
  u2'
  %  The next statement slows down the output.  You can increase the
  %  pause time by increasing the (0.1).
  
end