%%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

cfl=0.8;
cfl=input('CFL=   ');
%The number of mesh points and so on
I=400;
delx=1/(I-1);
i=1:I;
x=(i-1).*delx;


%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;

    
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(1:I-1)=1;            a1(I)=(1-cfl);
u11=uo';

left_matrix1=diag(a1)+diag(b1,-1)+diag(c1,1);
left_matrix_inv1=left_matrix1^-1;
%Second Derivative
u20=uo';
u21=u20;
u22=u20;
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_matrix2=diag(a2)+diag(b2,-1)+diag(c2,1)+diag(d2,-2);
left_matrix_inv2=left_matrix2^-1;

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

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