function lcrun
%LCRUN   Liquid crystal BVP.
%        Solves the liquid crystal BVP for four different lambda values.

lambda = [2.4, 2.5, 3, 10];

solinit = bvpinit(linspace(-1,1,20),@lcinit);
sola = bvp4c(@lc,@lcbc,solinit,[],lambda(1));
solb = bvp4c(@lc,@lcbc,solinit,[],lambda(2));
solc = bvp4c(@lc,@lcbc,solinit,[],lambda(3));
sold = bvp4c(@lc,@lcbc,solinit,[],lambda(4));
plot(sola.x,sola.y(1,:),'-','LineWidth',4), hold on
plot(solb.x,solb.y(1,:),'--','LineWidth',2)
plot(solc.x,solc.y(1,:),'--','LineWidth',4)
plot(sold.x,sold.y(1,:),'--','LineWidth',6), hold off
legend([repmat('\lambda = ',4,1) num2str(lambda')])
xlabel('x','FontSize',16)
ylabel('\theta','Rotation',0,'FontSize',16)

% ----------------------------------------------------------------------
% Subfunctions.
% ----------------------------------------------------------------------
function yprime = lc(x,y,lambda)
%LC       ODE/BVP liquid crystal system.
%         YPRIME = LC(X,Y,LAMBDA) evaluates derivative.

yprime = [y(2); -lambda*sin(y(1))*cos(y(1))];

% ----------------------------------------------------------------------
function res = lcbc(ya,yb,lambda)
%LCBC     ODE/BVP liquid crystal boundary conditions.
%         RES = LCBC(YA,YB,LAMBDA) evaluates residual.

res = [ya(1); yb(1)];

% ----------------------------------------------------------------------
function yinit = lcinit(x)
%LCINIT   ODE/BVP liquid crystal initial guess.
%         YINIT = LCINIT(X) evaluates initial guess at X.

yinit = [sin(0.5*(x+1)*pi); 0.5*pi*cos(0.5*(x+1)*pi)];
