diff --git a/chunkie/guide/addpaths_loc.m b/chunkie/guide/addpaths_loc.m
index 64bc844..92107db 100644
--- a/chunkie/guide/addpaths_loc.m
+++ b/chunkie/guide/addpaths_loc.m
@@ -1,4 +1,6 @@
 function addpaths_loc()
 
 addpath('../');
-addpath(genpath('../FLAM'));
\ No newline at end of file
+addpath('../fmm2d/matlab/');
+addpath('../FLAM');
+addpath(genpath_ex('../FLAM'));
\ No newline at end of file
diff --git a/chunkie/guide/guide_simplebvps.m b/chunkie/guide/guide_simplebvps.m
index 5d02360..ed83bb3 100644
--- a/chunkie/guide/guide_simplebvps.m
+++ b/chunkie/guide/guide_simplebvps.m
@@ -18,8 +18,8 @@
 
 % define a boundary condition. because of the symmetries of the 
 % starfish, this function integrates to zero
-f = @(r) r(1,:).^2.*r(2,:);
-rhs = f(chnkr.r); rhs = rhs(:);
+pwfun = @(r) r(1,:).^2.*r(2,:);
+rhs = pwfun(chnkr.r); rhs = rhs(:);
 
 % get the kernel 
 kernsp = kernel('lap','sprime');
@@ -46,9 +46,138 @@
 % plot
 figure(1); clf
 h = pcolor(xx,yy,uu); set(h,'EdgeColor','none'); colorbar
+hold on; plot(chnkr,'k'); colormap(redblue)
 % END LAPLACE NEUMANN
 saveas(figure(1),"guide_simplebvps_laplaceneumann.png")
 
+%%
 % START BASIC SCATTERING
+% get a chunker discretization of a peanut-shaped domain
+modes = [1.25,-0.25,0,0.5];
+ctr = [0;0];
+chnkr = chunkerfunc(@(t) chnk.curves.bymode(t,modes,ctr));
 
+% the boundary condition is determined by an incident plane wave
+kwav = 3*[-1,-5];
+pwfun = @(r) exp(1i*kwav*r(:,:));
+rhs = -pwfun(chnkr.r); rhs = rhs(:);
+
+% define the CFIE kernel D - ik S
+zk = norm(kwav);
+coefs = [1,-1i*zk];
+kerncfie = kernel('h','c',zk,coefs);
+
+% get a matrix discretization of the boundary integral equation 
+sysmat = chunkermat(chnkr,kerncfie); 
+% add the identity term
+sysmat = sysmat + 0.5*eye(chnkr.npt);
+
+% solve the system 
+sigma = gmres(sysmat,rhs,[],1e-10,100);
+
+% grid for plotting solution (in exterior)
+x1 = linspace(-5,5,300);
+[xx,yy] = meshgrid(x1,x1);
+targs = [xx(:).'; yy(:).'];
+in = chunkerinterior(chnkr,{x1,x1});
+uu = nan(size(xx));
+
+% same kernel to evaluate as on boundary
+uu(~in) = chunkerkerneval(chnkr,kerncfie,sigma,targs(:,~in));
+uu(~in) = uu(~in) + pwfun(targs(:,~in)).';
+
+% plot
+figure(1); clf
+h = pcolor(xx,yy,real(uu)); set(h,'EdgeColor','none'); colorbar
+colormap(redblue); umax = max(abs(uu(:))); clim([-umax,umax]);
+hold on 
+plot(chnkr,'k')
 % END BASIC SCATTERING
+saveas(figure(1),"guide_simplebvps_basicscattering.png")
+
+%%
+% START STOKES VELOCITY
+% get a chunker discretization of a peanut-shaped domain
+modes = [2.5,0,0,1];
+ctr = [0;0];
+chnkrouter = chunkerfunc(@(t) chnk.curves.bymode(t,modes,ctr));
+
+% inner boundaries are circles (reverse them to get orientations right)
+chnkrcirc = chunkerfunc(@(t) chnk.curves.bymode(t,0.3,[0;0]));
+chnkrcirc = reverse(chnkrcirc);
+centers = [ [-2:2, -2:2]; [(0.7 + 0.25*(-1).^(-2:2)) , ...
+    (-0.7 + 0.25*(-1).^(-2:2))]];
+centers = centers + 0.1*randn(size(centers));
+
+% make a boundary out of the outer boundary and several shifted circles
+chnkrlist = [chnkrouter];
+for j = 1:size(centers,2)
+    chnkr1 = chnkrcirc;
+    chnkr1.r(:,:) = chnkr1.r(:,:) + centers(:,j);
+    chnkrlist = [chnkrlist chnkr1];
+end
+chnkr = merge(chnkrlist);
+
+%%
+
+% boundary condition specifies the velocity. two values per node
+wid = 0.3; % determines width of Gaussian
+f = @(r) [exp(-r(2,:).^2/(2*wid^2)); zeros(size(r(2,:)))];
+rhsout = f(chnkrouter.r(:,:)); rhsout = rhsout(:);
+rhs = [rhsout; zeros(2*10*chnkrcirc.npt,1)];
+
+% define the combined layer Stokes representation kernel
+mu = 1; % stokes viscosity parameter
+kerndvel = kernel('stok','dvel',mu);
+kernsvel = kernel('stok','svel',mu);
+
+% get a matrix discretization of the boundary integral equation 
+dmat = chunkermat(chnkr,kerndvel); 
+smat = chunkermat(chnkr,kernsvel);
+
+% add the identity term and nullspace correction
+c = -1;
+W = normonesmat(chnkr);
+sysmat = dmat + c*smat - 0.5*eye(2*chnkr.npt) + W;
+
+% solve the system 
+sigma = gmres(sysmat,rhs,[],1e-10,100);
+
+% grid for plotting solution (in exterior)
+x1 = linspace(-3.75,3.75,300);
+y1 = linspace(-2,2,300);
+[xx,yy] = meshgrid(x1,y1);
+targs = [xx(:).'; yy(:).'];
+in = chunkerinterior(chnkr,{x1,y1});
+uu = nan([2,size(xx)]);
+pres = nan(size(xx));
+
+% same kernel to evaluate as on boundary
+uu(:,in) = reshape(chunkerkerneval(chnkr,kerndvel,sigma,targs(:,in)),2,nnz(in));
+uu(:,in) = uu(:,in) + c*reshape(chunkerkerneval(chnkr,kernsvel,sigma,targs(:,in)),2,nnz(in));
+uu(:,in) = uu(:,in) + uconst;
+
+% can evaluate the associated pressure
+kerndpres = kernel('stok','dpres',mu);
+kernspres = kernel('stok','spres',mu);
+opts = []; opts.eps = 1e-3;
+pres(in) = chunkerkerneval(chnkr,kerndpres,sigma,targs(:,in),opts);
+pres(in) = pres(in) + c*chunkerkerneval(chnkr,kernspres,sigma,targs(:,in),opts);
+
+% plot
+figure(1); clf
+h = pcolor(xx,yy,pres); set(h,'EdgeColor','none'); colorbar
+colormap("parula");
+hold on
+plot(chnkr,'k','LineWidth',2); axis equal
+
+figure(1)
+u = reshape(uu(1,:,:),size(xx)); v = reshape(uu(2,:,:),size(xx));
+startt = linspace(pi-pi/12,pi+pi/12,20);
+startr = 0.99*chnk.curves.bymode(startt,modes,[0;0]);
+startx = startr(1,:); starty = startr(2,:);
+
+sl = streamline(xx,yy,u,v,startx,starty);
+set(sl,'LineWidth',2,'Color','w')
+% END STOKES VELOCITY PROBLEM
+saveas(figure(1),"guide_stokesvelocity.png")
diff --git a/chunkie/guide/guide_simplebvps_laplaceneumann.png b/chunkie/guide/guide_simplebvps_laplaceneumann.png
index a2a8afd..22e1bc1 100644
Binary files a/chunkie/guide/guide_simplebvps_laplaceneumann.png and b/chunkie/guide/guide_simplebvps_laplaceneumann.png differ
diff --git a/docs/guide/simplebvps.rst b/docs/guide/simplebvps.rst
index c5d70fc..56d6bfc 100644
--- a/docs/guide/simplebvps.rst
+++ b/docs/guide/simplebvps.rst
@@ -228,13 +228,86 @@ can be solved with very little coding:
 
 .. image:: ../../chunkie/guide/guide_simplebvps_laplaceneumann.png
    :width: 500px
-   :alt: combined chunker
+   :alt: solution of Laplace equation
    :align: center
 		
 
 A Helmholtz Scattering Problem
 -------------------------------
 
+In a scattering problem, an incident field $u^{\\textrm{inc}}$ is specified
+in the exterior of an object and a scattered field $u^{\\textrm{scat}}$ is
+determined so that the total field $u = u^{\\textrm{inc}}+u^{\\textrm{scat}}$
+satisfies the PDE and a given boundary condition. It is also required that 
+the scattered field radiates outward from the object. For a sound-soft object
+in acoustic scattering, the total field is zero on the object boundary.
+This results in the following boundary value problem for $u^{\\textrm{scat}}$
+in the exterior of the object:
+
+.. math::
+
+   \begin{align*}
+   (\Delta + k^2) u^{\textrm{scat}} &= 0 & \textrm{ in } \mathbb{R}^2 \setminus \Omega \; , \\
+   u^{\textrm{scat}} &= -u^{\textrm{inc}} & \textrm{ on } \Gamma \; , \\
+   \sqrt{|x|} \left( \frac{x}{|x|} \cdot \nabla u^{\textrm{scat}} - ik u^{\textrm{scat}} \right )
+   &\to 0 & \textrm{ as } |x|\to \infty \; ,
+   \end{align*}
+
+The Green function for the Helmholtz equation is
+
+.. math::
+
+   G_k (x,y) = \frac{i}{4} H_0^{(1)}( k|x-y|) \; .
+
+Analogous to the above, this Green function can be used to define
+single and double layer potential operators
+
+.. math::
+
+   \begin{align*}
+   [S\sigma](x) &:= \int_\Gamma G_k(x,y) \sigma(y) ds(y)  \\
+   [D\sigma](x) &:= \int_\Gamma n(y)\cdot \nabla_y G_k(x,y) \sigma(y) ds(y) 
+   \end{align*}
+   
+A robust choice for the layer potential representation for this problem is
+a *combined field* layer potential, which is a linear combination
+of the single and double layer potentials:
+
+.. math::
+
+   u^{\textrm{scat}}(x) = [(D - ik S)\sigma](x) \; .
+
+Then, imposing the boundary conditions on this representation
+results in the equation
+
+.. math::
+
+   \begin{align*}
+   -u^{\textrm{inc}}(x_0) &= \lim_{x\in \mathbb{R}^2\setminus \Omega, x\to x_0} [(D-ik S)\sigma](x)  \\
+   &= \frac{1}{2} \sigma(x_0) + [ (\mathcal{D} -ik \mathcal{S})\sigma ](x_0)
+   \end{align*}
+
+where we have used the exterior jump condition for the double layer potential.
+As above, the integrals in the operators restricted to the boundary must sometimes
+be interpreted in the principal value or Hadamard finite-part sense.
+The above is a second kind integral equation and is invertible.
+
+
+As before, this is relatively straightforward to implement in :matlab:`chunkie`
+using the built-in kernels:
+
+.. include:: ../../chunkie/guide/guide_simplebvps.m
+   :literal:
+   :code: matlab
+   :start-after: % START BASIC SCATTERING
+   :end-before: % END BASIC SCATTERING
+
+
+.. image:: ../../chunkie/guide/guide_simplebvps_basicscattering.png
+   :width: 500px
+   :alt: acoustic scattering around a peanut shape
+   :align: center
+