clear all;clc;

global X1 elcon;

%*********** Geometry ***********************************
QL = [1];

%*********** Physical Parameters ***********************
Xleft = 0;
Xright = 40;
nnodes = 41;
X1 = linspace(Xleft,Xright,nnodes);
u = ones(nnodes,1);
delta_x = (X1(length(X1))-X1(1))/(nnodes-1);
qleft = 0;
nodeL = ceil(length(X1)/4);
L = X1(nodeL)-X1(1);
XL = X1(1:nodeL);
q0 = zeros(length(X1),1);
q0(1:nodeL) = (1-cos(2*pi*XL/L))/2;

%*********** Time Integration Parameters ***********************
courant = 0.5;
delta_t = courant*delta_x/u(1);
theta = 0.5;
start = 1;
inc = 1;
stop = 61;
[Xnodes Ynodes] = size(X1);
q = q0;
q_last = q;
%unsteady time integration loop.
time_count = 1;
for timeLoop = start:inc:stop;
q = q_last;

%** Apply Boundary Conditions to Initial Condition for Set 0 ********************
for ii=1:1:length(QL)
q(QL(ii)) = +qleft;
end

%********************* Equations for Set 0*****************
q0=asres1D(delta_t,[],u,0,'A3001L',q);

%Apply Dirichlet boundary condition to boundary QL.
for ii=1:1:length(QL)
q0(QL(ii)) = 0;
end

%***************** Jacobians for Set 0**********
%Equation q by q
q_jac0 = +asjac1D(1,[],[],1,'A200L',[])...
+asjac1D(theta.*delta_t,[],u,0,'A3001L',[]);

%Apply Dirichlet boundary condition to boundary QL.
for ii=1:1:length(QL)
q_jac0(QL(ii),:) = 0;
q_jac0(:,QL(ii)) = 0;
q_jac0(QL(ii),QL(ii)) = 1;
end

%***************** Solver ********************************************
dq = q_jac0\(-q0);

%************* Solve Iterate *****************************************
q = q + dq;
q_store(:,time_count) = q;
figure(1)
hold on
q_last = q;
time_count = time_count + 1;
if time_count==2    
  plot(X1,q,'k')
elseif time_count==21
  plot(X1,q,'g')
elseif time_count==41
  plot(X1,q,'r')
elseif time_count==61
  plot(X1,q,'c')
  %elseif time_count==41
  %plot(X1,q,'m')
  %elseif time_count==51
  %plot(X1,q,'y')
  xlabel('X1','FontSize',14),ylabel('q','FontSize',14)
  title(['Courant=' num2str(courant) ', Theta=' num2str(theta)] , 'FontSize',14)
  legend('t=1','t=20','t=40','t=60','Location',([0.62 .715 0.161 0.194]))
else
end 
end %end time integration loop
