%%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 implicit time marching methods.
%
%The initial condition is specified and the wave speed is assumed to be
%positive
%
%For nonlinear case
%CFL=u*delt/delx;	-->delt/delx=CFL/u

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

delxarr=repmat(delx,[1 I]);
cfl=1.8;
pause
%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;

nsteps=10;
%METHOD 1-----------------------Non-linear backwards difference

u=uo;
for i=1:nsteps
    f=u.^2./2;                  %u du/dx=.5* d(u^2)/dx
    delt=cfl.*delxarr./u;
    u(2:I)=u(2:I)-(delt(2:I)./delx).*(u(2:I)).*(u(2:I)-u(1:I-1));
    %u(2:I)=u(2:I)-(delt(2:I)./delx).*(f(2:I)-f(1:I-1));
    plot(x,u);
    pause;
    set(gca,'ylim',[.25,1.25]);
end


u=uo;
for i=1:nsteps
   delt=cfl.*delxarr./u
   a=zeros(1,I);                                    a(1:I-1)=1;    a(I)=1+(delt(I)/delx)*u(I);
   b(1:I-2)=-.5*(delt(2:I-1)./delx).*u(1:I-2);                     b(I-1)=-(delt(I)/delx)*u(I-1);
   c(2:I-1)= .5*(delt(2:I-1)./delx).*u(3:I);                       c(1)=0;

   [L,U,alph,gamm]=LU_decomposition(a,b,c,I);
   f(1)=uo(1);              f(I)=u(I)+.5*(delt(I)./delx).*(u(I)^2-u(I-1)^2);
   f(2:I-1)=u(2:I-1)+.25*(delt(2:I-1)/delx).*(u(3:I).^2-u(1:I-2).^2);
   u=LU_solve(alph,b,gamm,f,I)';
   plot(x,u);
   set(gca,'ylim',[0.25,1.25]);
   pause;
end