clear all;clc;

global X1 elcon;
nnodes=21;

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

% Code to generate non-uniform mesh using Newton-Raphson 
% Iterative method ****
iter=0;
itermax=20;
eps=1e-16;
a=0.95;
r=0.7; % progression ratio
S=0;
dx=abs(S-1);
while dx > eps & iter < itermax
S=0;
den=0;
M=nnodes-1;
for ii=1:M
    n=ii+1;
    x(n)=a*r.^(n-1);
    S=x(n)+S;
    den=den+r^(ii-1);
end

num=S-1;
anew=a-(num)/den;
dx=abs(anew-a);
a=anew;
iter = iter + 1;
end
x;
X1=x;
for i=1:(M+1)
    X1(i)=X1(i);
end
for i=2:(M+1)
    X1(i)=X1(i)+X1(i-1);
end
% Non-uniform mesh generation ends

%*********** Physical Parameters ***********************
% Uniform mesh 
%left = 0;
%right = 1.0;
%nnodes = 17;
%X1 = linspace(left,right,nnodes);
one_Pe = 0.001;
Pe=1/one_Pe;
u = ones(nnodes,1);
Qleft = 0;
Qright = 1;
[Xnodes Ynodes] = size(X1);
q = ones(Xnodes.*Ynodes,1)*0;

%********************* 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(1,[],u,0,'A3001L',q)...
+asres1D(one_Pe,[],[],-1,'A211L',q);

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

%Apply Neumann boundary condition to boundary QR.
q0(nnodes) = q0(nnodes) -Qright;

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

%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
figure(2)

plot(X1,q,'md-')
%title('Scalar Field Transport, Non-Uniform Mesh, Pe=1000','FontSize',16)
title({'Scalar Field Transport, Non-Uniform Mesh, Pe=' num2str(Pe)},'FontSize',16)
%xlabel('X1','FontSize',16),ylabel('q','FontSize',16)
xlabel('X1','FontSize',16),ylabel('q','FontSize',16)
%Axis(0,1,0,1)