%%Lance Simms                                               2/13/06
%%AA214B-MacCormack
%
%This Script will attempt to solve the model hyperbolic equation 
%
%       du    du
%       -- +c -- = 0
%       dt    dx
%
%using several explicit time marching methods.
%
%The initial condition is specified and the wave speed is assumed to be
%positive

%The number of mesh points and so on
I=41;
delx=2/(I-1);
i=1:I;
x=(i-1).*delx;
cfl=.9;

%Initial condition has first half of wave at 1 and second half at 1/2
uo=zeros(1,I);
uo(i(find(i < I/4+1)))=1;
uo(i(find(i > I/4+1)))=1/2;

%1st Method------Bacward Explicit
u=uo;
nsteps=10;

if im_ex==1
    
for i=1:nsteps  
    u(2:I)=u(2:I)-cfl*(u(2:I)-u(1:I-1));
    %plot(x,u);%pause;
    
end
    subplot(4,2,1);
    plot(x,u);
    set(gca,'units','normalized');
    text(.6,.8,'Explicit Backwards','units','normalized');

    u=uo;
%2nd Method------Forward Explicit
for i=1:nsteps  
    u(2:I-1)=u(2:I-1)-cfl*(u(3:I)-u(2:I-1));
    %plot(x,u)%;pause;
    
end

    subplot(4,2,2);
    plot(x,u);
    text(.6,.8,'Explicit Forward','units','normalized');

    
u=uo
for i=1:nsteps  
    u(2:I-1)=u(2:I-1)-cfl*(u(3:I)-u(1:I-2))/2;
    u(I)=u(I)-cfl*(u(I)-u(I-1));
    %plot(x,u)%pause; 
end

    subplot(4,2,3);
    plot(x,u);
    text(.6,.8,'Explicit Central','units','normalized');

%6th Method--------Lax
u=uo
for i=1:nsteps  
    u(2:I-1)=(u(3:I)+u(1:I-2))/2-cfl*(u(3:I)-u(1:I-2))/2;
    u(I)=u(I)-cfl*(u(I)-u(I-1));
    %plot(x,u);pause;
end

    subplot(4,2,4);
    plot(x,u);
    text(.6,.8,'Lax','units','normalized');

%7th Method--------Lax-Wendroff
u=uo
for i=1:nsteps 
    u(2:I-1)=u(2:I-1)-cfl*(u(3:I)-u(1:I-2))/2+...
            .5*cfl^2*(u(3:I)-2*u(2:I-1)+u(1:I-2));

    u(I)=u(I)-cfl*(u(I)-u(I-1));
    %plot(x,u); pause;
end

    subplot(4,2,5);
    plot(x,u);
    text(.6,.8,'Lax-Wendroff','units','normalized');

%8th Method--------MacCormack
u=uo
ubar=uo
for i=1:nsteps 
    ubar(2:I-1)=u(2:I-1)-cfl*(u(3:I)-u(2:I-1));
    ubar(I)=ubar(I)-cfl*(ubar(I)-ubar(I-1));
    u(2:I-1)=.5*(u(2:I-1)+ubar(2:I-1)-...
             cfl*(ubar(2:I-1)-ubar(1:I-2)));
    u(I)=u(I)-cfl*(u(I)-u(I-1));
    %plot(x,u); pause;
end

    subplot(4,2,6);
    plot(x,u);
    text(.6,.8,'MacCormack','units','normalized');

%9th Method--------Jameson
u=uo
for i=1:nsteps 
    uk=u;
    for k=1:4
        alphak=1/(5-k);
        uk(2:I-1)=u(2:I-1)-alphak*.5*cfl*(uk(3:I)-uk(1:I-2));
        uk(I)=uk(I)-cfl*(uk(I)-uk(I-1));
    end
    u=uk
    %plot(x,u); pause;
end

    subplot(4,2,7);
    plot(x,u);
    text(.6,.8,'Jameson','units','normalized');
    text(1.4,.8,'X-axis: x','units','normalized');
    text(1.4,.65,'Y-axis: u (Displacement)','units','normalized');
    text(1.4,.4,sprintf('CFL= %6.2f', cfl),'units','normalized')

end

if im_ex==2

%Central Implicit Method
   u=uo;
   alpha=1;
   a=zeros(1,I);   a(1:I-1)=1;      a(I)=1+cfl;
   b(1:I-2)=-.5*alpha*cfl;          b(I-1)=-cfl;
   c(2:I-1)=-.5*alpha*cfl;          c(1)=0;
   [L,U,alph,gamm]=LU_decomposition(a,b,c,I);
   
   for i=1:nsteps
       f(1)=uo(1);                              f(I)=u(I);
       f(2:I-1)=u(2:I-1)-(1-alpha)*cfl.*(u(3:I)-u(1:I-2));
       u