Function files bratuf.m, bratujac.m and bratuinit.m

function [c,a,f,b]=bratuf(p,u,lam) 
% pde for Bratu 
u=pdeintrp(p.points,p.tria,u); c=1; a=0;f=-10*(u-lam*exp(u)); b=0;

function [c,fu,flam,b]=bratujac(p,u,lam) 
% jacobian for Bratu 
u=pdeintrp(p.points,p.tria,u); c=1; 
fu=-10*(1-lam*exp(u)); flam=10*exp(u); b=0;

function p=bratuinit(p) 
% init-routine 
p=stanparam(p); % set standard parameters, usually call this first
p.neq=1; p.f=@bratuf; p.jac=@bratujac; % the coefficients and the jacobian
[p.geo,bc]=recnbc1(0.5,0.5); % geometry and boundary conditions 
p.bcf=@(p,u,lam) bc; % define dummy BC functions (as BC are fixed)
nx=20;p=stanmesh(p,nx,nx); % rectangular mesh 
p=setbmesh(p); % set current mesh as 'base mesh' for adaption 
p.maxt=2*p.nt; % usual goal for mesh-adaption 
pre=sprintf('%s',inputname(1)); p=setfn(p,pre); % set filenames  
p.tau=zeros(p.neq*p.np+1,1); p.tau(1)=1; % init first tangent 
p.xi=1/p.np; % and set xi to its standard value 
p.pstyle=1; p.cm=hot; % mesh-plot as standard, colormap for 3D surf plot
p.nsteps=50; p.dsmax=1; % number of steps and maximal stepsize 
p.dlammax=0.02; % maximal stepsize in lambda; 
p.lammin=0.05;  % minimal lambda, used as stopping criterion
p.lam=0.2; % starting point in lam, and in u: 
p.u=0.1*ones(p.np,1);p.ds=0.05;  % u on lower homogen. branch 
% next (and final) steps not necessary but good practice to assess 
% quality of initial guess and provide first point on branch
res=norm(resi(p,p.u,p.lam), p.normsw); % calculate first residual 
fprintf('initial res=%g\n',res); % and print it 
[p.u,p.res,p.iter]=nlooppde(p,p.u,p.lam); % correct initial guess 
fprintf('first res=%g with %i iterations\n',p.res,p.iter); % and inform user 
plotsol(p,1,1,p.pstyle);

Script bratucmds.m

% script 
close all; clear all; path('../../pde2p',path); path('../../aux',path);
p=[]; p=bratuinit(p); p=cont(p); % homogeneous branch 
%% 2 branch switchings and successive continuation
q=swibra('p','bp1','q'); q.lammin=0.09;q.nsteps=21;q=cont(q); 
r=swibra('q','bp1','r'); r.lammin=0.05;r.nsteps=20;r=cont(r); 
%% postprocessing, first plot BifDiagram 
figure(3);clf(3);
plotbra(p,3,2,'ms',12,'lwst',5,'lwun',2,'fs',16,'cl','k');
plotbra(q,3,2,'ms',12,'lwst',6,'lwun',2,'lab',[20],'fs',16,'cl','b');
plotbra(r,3,2,'ms',12,'lwst',6,'lwun',2,'lab',[20],'fs',16,'cl','r');
axis([0.05 0.37 0.5 3.5]);xlabel('\lambda');ylabel('||u||_2');
%% plot solutions 
plotsolf('q','p20',7,1,1); % mesh plot 
figure(7); view(-30,50); xlabel('x');ylabel('y');
plotsolf('r','p20',8,1,3); % rendered 3D plot 
figure(8); view(-30,50); xlabel('x');ylabel('y');