function GradKansa(M,N,ep)
% Exact solution
u = @(x,ep) (1+exp(-1/ep)-exp(-x/ep)-exp((x-1)/ep));
% xx der of Exact solution
uxx = @(x,ep) -(exp(-x/ep)/ep^2 + exp((x - 1)/ep)/ep^2);
% RHS
fu = @(x,ep) -ep^2*uxx(x,ep)-u(x,ep);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
rbf = @(r,c) sqrt((r*c).^2+1);
rbfxx = @(r,c)  c^2./((r*c).^2+1).^(3/2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%source pts
Nt = 51; % # of interior test points
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Generate collocation points
xp(1:M+1,1)=(0:M)./M;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Generate centres 
xs(1:N+1,1)=(0:N)./N;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Interior test points
xt(1:Nt+1,1)=(0:Nt)./Nt;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Formulate the matrix [ HL B1; H  B2; F  0] with augmented polynomial
DM1(1:M+1,1:N+1)=xp(1:M+1,1).*ones(1,N+1)-ones(M+1,1).*(xs(1:N+1,1))';
DM2(1:Nt+1,1:N+1)=xt(1:Nt+1,1).*ones(1,N+1)-ones(Nt+1,1).*(xs(1:N+1,1))';
%% RHS
rhs1 = zeros(M+1,1);
rhs1(1,1)=u(xp(1,1),ep);
rhs1(2:M,1)=fu(xp(2:M,1),ep);
rhs1(M+1,1)=u(xp(M+1,1),ep);
%% Global matrix
for iter=1:100
c=0.1*iter;    
G1=rbf(DM1,c);
G2=rbfxx(DM1,c);
GL(1,1:N+1)=G1(1,1:N+1);
GL(2:M,1:N+1)= -ep^2*G2(2:M,1:N+1)-G1(2:M,1:N+1);
GL(M+1,1:N+1)=G1(M+1,1:N+1);
%GL
%%% Solution of system
sol=GL\rhs1;
% compute the approximate solution
B = rbf(DM2,c);
uh = B*sol;    % approximate solution
ut = u(xt,ep);  % exact solution
%uh-ut
error = norm(uh-ut,inf)/norm(ut,inf);  %relative error
rmse=norm(uh-ut,2)/sqrt(length(xt));   %RMSE error
maxer=norm(uh-ut,inf);  % maximum error
fprintf('epsilon=%5.2f, M=%4d, N=%4d, c =%5.3f,  rel_er=%10.3e, RMSE=%10.3e\n',ep,M,N,c,error,rmse) 
xxx(iter)=c;yyy(iter)=error;zzz(iter)=rmse;
end
semilogy(xxx,yyy,'b-',xxx,zzz,'r-.','LineWidth',2)
%xlim([0 1]);%set(gca,'xtick',[ 0 0.05  0.1]);
%xlim([0  50]);set(gca,'xtick',[ 0 10  20  30  40  50]);
title('Errors vs c')
legend('Rel Error','RMSE','Location','NorthEast')
ylabel('Error'),xlabel('c');
%fprintf('epsilon=%5.2f, rel_er=%10.3e, max_er=%10.3e, RMSE=%10.3e\n',ep,error,maxer,rmse);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% End of the main program