How to solve 2nd order BVP using Chebyshev Spectral Method?

Hi, I am tring to solve the following BVP using Chebshev Spectral Method:

with boundary conditions:

I am following the code written by Lloyd N. Trefethen in his book "SPECTRAL METHODS IN MATLAB".

The code is given by for Linear BVP

% p33.m - solve linear BVP u_xx = exp(4x), u'(-1)=u(1)=0

N = 16;
[D,x] = cheb(N); D2 = D^2;              
D2(N+1,:) = D(N+1,:);            % Neumann condition at x = -1
D2 = D2(2:N+1,2:N+1);                
f = exp(4*x(2:N));           
u = D2\[f;0];                       
u = [0;u];              
clf, subplot('position',[.1 .4 .8 .5])
plot(x,u,'.','markersize',16)
axis([-1 1 -4 0])
xx = -1:.01:1;
uu = polyval(polyfit(x,u,N),xx);   
line(xx,uu)
grid on
exact = (exp(4*xx) - 4*exp(-4)*(xx-1) - exp(4))/16; 
title(['max err = ' num2str(norm(uu-exact,inf))],'fontsize',12)

where cheb was defined as:

% CHEB  compute D = differentiation matrix, x = Chebyshev grid

function [D,x] = cheb(N)
if N==0, D=0; x=1; return, end
x = cos(pi*(0:N)/N)'; 
c = [2; ones(N-1,1); 2].*(-1).^(0:N)';
X = repmat(x,1,N+1);
dX = X-X';                  
D  = (c*(1./c)')./(dX+(eye(N+1)));      % off-diagonal entries
D  = D - diag(sum(D'));                 % diagonal entries

I want to modify the the code according to the BVP defined earlier. Please help me to sort out this problem.

Thanks

[Torsten]

You are forced to use this special solution method ? Or are you also free to use bvp4c or bvp5c ?

[Muhammad Usman]

I want to replicate the result generated in an research article with the method I mentioned in my question and expected result is as follows:

Answer 1

[Torsten]

alpha = 0.1;
xmesh = linspace(0,1,50);
solinit = bvpinit(xmesh, [0 0]);
bvpfcn = @(x,y)[y(2);0.75/(1+6*alpha*y(2)^2)];
bcfcn = @(ya,yb)[ya(2);yb(1)];
sol = bvp4c(bvpfcn, bcfcn, solinit);
plot(sol.x, sol.y, '-o')
grid on

[Muhammad Usman]

This is not the expected result, I showed the expected result above in the comment

[Torsten]

This curve does not satisfy your boundary conditions.

Further, it seems to have a vertical tangent at its x-crossing - thus u' = Inf there.