%%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=100;
delx=1/(I-1);
i=1:I;
x=(i-1).*delx;
cfl=-1;

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

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

cfl=8;


    
u0=uo;
u1=uo;
u2=uo;

%%Explicit Backwards
for j=1:nsteps  
    
    u2(2:I-1)=2*u1(2:I-1)-u0(2:I-1)-cfl^2.*(u1(3:I)-2*u1(2:I-1)+u1(1:I-2));
    u2(I)=2*u1(I)-u0(I)+cfl^2.*(u1(I)-2*u1(I-1)+u1(I-2));
    u1=u2;
    u0=u1;
    %plot(x,u2);pause; 
end

u0=uo;
u1=uo;
u2=uo;

%%Explicit Backwards
for j=1:nsteps  
    
    u2(3:I)=2*u1(3:I)-u0(3:I)+cfl^2.*(u1(3:I)-2*u1(2:I-1)+u1(1:I-2));
    u2(2)=2*u1(2)-u0(2)+cfl^2.*(u1(3)-2*u1(2)+u1(1));
    u1=u2;
    u0=u1; 
    %plot(x,u2);pause; 
    
end

%%Implicit Central
initial_condition=1
if initial_condition==1
    uo(:)=0
end

%First Derivative
b1(1:I-2)=-cfl;         b1(I-1)=cfl;
c1(2:I-1)=-cfl;         c1(1)=0;
a1(2:I-1)=
%Second Derivative
u20=uo';
u21=u0;
u22=u0;
b2(1:I-2)=-cfl^2;        b2(I-1)=2*cfl^2;
c2(2:I-1)=-cfl^2;        c2(1)=(0);
d2(1:I-3)=0;             d2(I-2)=-cfl^2;
a2(2:I-1)=(1+2*cfl^2);   a2(1)=1;         a2(I)=(1-cfl^2);

left_matrix=diag(a2)+diag(b2,-1)+diag(c2,1)+diag(d2,-2);
left_matrix_inv=left_matrix^-1;

plot(uo);
set(gcf,'doublebuffer','on');
nsteps=800

for j=1: nsteps
    f2=2*u21-u20;
    f2(1)=.5*sin(2*pi*j/80);
    u22=left_matrix_inv*f2;
    u21=u22;
    u20=u21; 
    plot(x,u22);
    set(gca,'ylim',[-2,2]);
    pause;
    %Form the left hand matrix a main diagonal, b below, c above
end