function [x,y,N]=BLGrid(bl, N_bl)
%x,y grid
%N total number of nodes
% bl boundary layer thickness
%N_bl the number of nodes in the boundary layer region


x1=zeros(N_bl,1);
hunif=(bl^(1/2))/(N_bl-1);
for i=1:N_bl
    x1(i)=((i-1)*hunif)^2;
end
hmin=x1(2)-x1(1);
hmax=(x1(N_bl)-x1(N_bl-1));
Nm=floor(((1-bl-(bl))/hmax))-1;
N=Nm+2*N_bl-1;
h=(1-bl-(bl))/Nm;
x2=zeros(Nm-1,1);
for i=1:Nm-1
    x2(i)=bl+i*h;
end

x3=zeros(N_bl,1);
hunif=(1^(1/2)-(1-bl)^(1/2))/(N_bl-1);
x3(N_bl)=1;
for i=N_bl-1:-1:1
    x3(i)= 1-x1(N_bl-i+1);
end
xunif=0:1/(N-1):1;

x=[x1;x2;x3];
y=x;
[N,M]=size(x);

% figure(1)
% plot(x,xunif)

% figure(3)
% hold on
% for i=1:N
% plot(x,ones(N,1)*y(i));
% end
% for i=1:N
% plot(ones(N,1)*x(i),y);
% end
% axis equal
% axis off