diff --git a/chunkie/+chnk/+axissymhelm2d/axissymhelm2d.m b/chunkie/+chnk/+axissymhelm2d/axissymhelm2d.m
new file mode 100644
index 00000000..05b76244
--- /dev/null
+++ b/chunkie/+chnk/+axissymhelm2d/axissymhelm2d.m
@@ -0,0 +1,108 @@
+function obj = axissymhelm2d(type, zk, coefs)
+%KERNEL.AXISSYMHELM2D   Construct the axissymmetric Helmholtz kernel.
+%   KERNEL.AXISSYMHELM2D('s', ZK) or KERNEL.AXISSYMHELM2D('single', ZK) 
+%   constructs the axissymmetric single-layer Helmholtz kernel with 
+%   wavenumber ZK.
+%
+%   KERNEL.AXISSYMHELM2D('d', ZK) or KERNEL.AXISSYMHELM2D('double', ZK) 
+%   constructs the axissymmetric double-layer Helmholtz kernel with 
+%   wavenumber ZK.
+%
+%   KERNEL.AXISSYMHELM2D('sp', ZK) or KERNEL.AXISSYMHELM2D('sprime', ZK) 
+%   constructs the normal derivative of the axissymmetric single-layer 
+%   Helmholtz kernel with wavenumber ZK.
+%
+%   KERNEL.AXISSYMHELM2D('c', ZK, COEFS) or 
+%   KERNEL.AXISSYMHELM2D('combined', ZK, COEFS)
+%   constructs the combined-layer axissymmetric Helmholtz kernel with 
+%   wavenumber ZK and parameter COEFS, 
+%   i.e., COEFS(1)*KERNEL.AXISSYMHELM2D('d', ZK) + 
+%         COEFS(2)*KERNEL.AXISSYMHELM2D('s', ZK).
+%
+% NOTES: The axissymetric kernels are currently supported only for purely
+% real or purely imaginary wave numbers
+%
+% See also CHNK.AXISSYMHELM2D.KERN.
+
+if ( nargin < 1 )
+    error('Missing Helmholtz kernel type.');
+end
+
+if ( nargin < 2 )
+    error('Missing Helmholtz wavenumber.');
+end
+
+
+zr = real(zk); zi = imag(zk);
+if abs(zr)*abs(zi) > eps
+    error('Only purely real or purely imaginary wavenumbers supported');
+end
+
+obj = kernel();
+obj.name = 'axissymhelmholtz';
+obj.params.zk = zk;
+obj.opdims = [1 1];
+
+switch lower(type)
+
+    case {'s', 'single'}
+        obj.type = 's';
+        obj.eval = @(s,t) chnk.axissymhelm2d.kern(zk, s, t, [0,0], 's');
+        obj.shifted_eval = @(s,t,o) chnk.axissymhelm2d.kern(zk, s, t, o, 's');
+        obj.fmm = [];
+        obj.sing = 'log';
+
+    case {'d', 'double'}
+        obj.type = 'd';
+        obj.eval = @(s,t) chnk.axissymhelm2d.kern(zk, s, t, [0,0], 'd');
+        obj.shifted_eval = @(s,t,o) chnk.axissymhelm2d.kern(zk, s, t, o, 'd');
+        obj.fmm = [];
+        obj.sing = 'log';
+
+    case {'sp', 'sprime'}
+        obj.type = 'sp';
+        obj.eval = @(s,t) chnk.axissymhelm2d.kern(zk, s, t, [0,0], 'sprime');
+        obj.shifted_eval = @(s,t,o) chnk.axissymhelm2d.kern(zk, s, t, o, 'sprime');
+        obj.fmm = [];
+        obj.sing = 'log';
+
+    case {'c', 'combined'}
+        if ( nargin < 3 )
+            warning('Missing combined layer coefficients. Defaulting to [1,1i].');
+            coefs = [1,1i];
+        end
+        obj.type = 'c';
+        obj.params.coefs = coefs;
+        obj.eval = @(s,t) chnk.axissymhelm2d.kern(zk, s, t, [0,0], 'c', coefs);
+        obj.shifted_eval = @(s,t,o) chnk.axissymhelm2d.kern(zk, s, t, o, 'c', coefs);
+        obj.fmm = [];
+        obj.sing = 'log';
+     case {'neu_rpcomb'}
+        obj.type = 'neu_rpcomb';
+        if ( nargin < 3 )
+            warning('Missing coefficient of 1i*D. Defaulting to 1.');
+            coefs = 1;
+        end
+        obj.eval = @(s,t) chnk.axissymhelm2d.kern(zk, s, t, [0,0], 'neu_rpcomb', coefs);
+        obj.shifted_eval = @(s,t,o) chnk.axissymhelm2d.kern(zk, s, t, o, 'neu_rpcomb', coefs);
+        obj.fmm = [];
+        obj.sing = 'log';
+        obj.opdims = [3 3];
+        obj.params.c1 = -1.0/(0.5 + 0.25*1i*coefs);
+        obj.params.c2 = -1i*coefs/(0.5 + 0.25*1i*coefs);
+
+    case {'dp_diff_d0p'}
+        obj.type = 'dp_diff_d0p';
+        obj.eval = @(s,t) chnk.axissymhelm2d.kern(zk, s, t, [0,0], 'dprime') - ...
+                            chnk.axissymlap2d.kern(s, t, [0,0], 'dprime', n);
+        obj.shifted_eval = @(s,t,o) chnk.axissymhelm2d.kern(zk, s, t, o, 'sprime') - ...
+                            chnk.axissymlap2d.kern(s, t, o, 'dprime', n);
+        obj.fmm = [];
+        obj.sing = 'log';
+
+    otherwise
+        error('Unknown axissym Helmholtz kernel type ''%s''.', type);
+
+end
+
+end
diff --git a/chunkie/+chnk/+axissymhelm2d/kern.m b/chunkie/+chnk/+axissymhelm2d/kern.m
index 213fa488..7463b109 100644
--- a/chunkie/+chnk/+axissymhelm2d/kern.m
+++ b/chunkie/+chnk/+axissymhelm2d/kern.m
@@ -299,6 +299,56 @@
   submat(1:3:end, 3:3:end) = c2*spikmat;
   submat(2:3:end, 1:3:end) = -sikmat;
   submat(3:3:end, 1:3:end) = -spikmat;
+
+case{'spp','sprimeprime'}
+  targnorm = targinfo.n(:,:);
+  [~,~,hess] = chnk.axissymhelm2d.green(zk, src, targ, origin);
+  nxtarg = repmat((targnorm(1,:)).',1,ns);
+  nytarg = repmat((targnorm(2,:)).',1,ns);
+  submat = hess(:,:,1).*nxtarg.*nxtarg + 2*hess(:,:,5).*nytarg.*nxtarg ...
+                                        + hess(:,:,3).*nytarg.*nytarg;
+
+  fker = @(x, s, t, rnt) fsprimeprime(x, zk, s, t, rnt, origin);
+    for j=1:ns
+        for i=1:nt
+            rt = targ(1,i) + origin(1);
+            dr = (src(1,j) - targ(1,i));
+            dz = (src(2,j) - targ(2,i));
+            r0   = sqrt(rt^2+(rt+dr)^2+dz^2);
+            alph = (dr^2+dz^2)/r0^2;
+            if alph > 2e-4 && alph < 0.2
+                [x0, w0] = get_grid(zk, rt, dr, dz);
+                fvals = fker(x0, src(:, j), targ(:,i), targnorm(:,i));
+                submat(i,j) = 2*w0.'*fvals;
+            end
+        end
+    end
+
+case {'q','quad','quadruple','quadrupole'} % q := spp'
+  srcnorm = srcinfo.n(:,:);
+  [~,~,hess] = chnk.axissymhelm2d.green(zk, src, targ, origin);
+  nxsrc = repmat(srcnorm(1,:),nt,1);
+  nysrc = repmat(srcnorm(2,:),nt,1);
+  submat = hess(:,:,2).*nxsrc.*nxsrc - 2*hess(:,:,6).*nysrc.*nxsrc ... 
+                                      + hess(:,:,3).*nysrc.*nysrc;
+
+  fker = @(x, s, t, rns) fqlp(x, zk, s, t, rns, origin);
+    for j=1:ns
+        for i=1:nt
+            rt = targ(1,i) + origin(1);
+            dr = (src(1,j) - targ(1,i));
+            dz = (src(2,j) - targ(2,i));
+            r0   = sqrt(rt^2+(rt+dr)^2+dz^2);
+            alph = (dr^2+dz^2)/r0^2;
+            if alph > 2e-4 && alph < 0.2
+                [x0, w0] = get_grid(zk, rt, dr, dz);
+                fvals = fker(x0, src(:, j), targ(:,i), srcnorm(:,i));
+                submat(i,j) = 2*w0.'*fvals;
+            end
+        end
+    end
+
+  
 otherwise
     error('Unknown axissymmetric Helmholtz kernel type ''%s''.', type); 
 end
@@ -357,6 +407,36 @@
            rndt.*rnds.*(-zk^2.*r.^2 - 3*1j*zk.*r + 3)./r.^5).*exp(1j*zk*r)/4/pi.*(rs + o(1));
 end
 
+function f = fqlp (x, zk, s, t, rns, o)
+    rs = s(1); zs = s(2);
+    rt = t(1); zt = t(2);
+
+    sxhalf = sin(x/2);
+    sxhalf2 = sxhalf.*sxhalf;
+    cx = 1-2*sxhalf2;
+    
+    rnds = ((rt + o(1)).*cx  - (rs + o(1))).*rns(1) + (zt - zs).*rns(2);
+    
+    r = sqrt((rs-rt).^2 + (zs-zt).^2 + 4*(rs+o(1)).*(rt+o(1)).*sxhalf2);
+    f = ((1j*zk.*r-1)./r.^3 + ...
+           rnds.^2.*(-zk.*r.^2 - 3*1j*zk.*r + 3)./r.^5).*exp(1j*zk*r)/4/pi.*(rs+o(1));
+end
+
+function f = fsprimeprime (x, zk, s, t, rnt, o)
+    rs = s(1); zs = s(2);
+    rt = t(1); zt = t(2);
+
+    sxhalf = sin(x/2);
+    sxhalf2 = sxhalf.*sxhalf;
+    cx = 1-2*sxhalf2;
+    
+    rndt = ((rt + o(1)) - (rs + o(1)).*cx).*rnt(1) + (zt - zs).*rnt(2);
+
+    r = sqrt((rs-rt).^2 + (zs-zt).^2 + 4*(rs+o(1)).*(rt+o(1)).*sxhalf2);
+    f = ((1j*zk.*r-1)./r.^3 + ...
+           rndt.^2.*(-zk.*r.^2 - 3*1j*zk.*r + 3)./r.^5).*exp(1j*zk*r)/4/pi.*(rs+o(1));
+end
+
 function [fkp, fik, fikp, fkdiff] = get_neu_kers(zk, cx, sxhalf2, s, t, rns, rnt, o)
 
     rs = s(1); zs = s(2);
diff --git a/chunkie/+chnk/+axissymlap2d/axissymlap2d.m b/chunkie/+chnk/+axissymlap2d/axissymlap2d.m
new file mode 100644
index 00000000..3fef1f36
--- /dev/null
+++ b/chunkie/+chnk/+axissymlap2d/axissymlap2d.m
@@ -0,0 +1,99 @@
+function obj = axissymlap2d(type, n)
+%KERNEL.AXISSYMHELM2D   Construct the axissymmetric Helmholtz kernel.
+%   KERNEL.AXISSYMHELM2D('s', ZK) or KERNEL.AXISSYMHELM2D('single', ZK) 
+%   constructs the axissymmetric single-layer Helmholtz kernel with 
+%   wavenumber ZK.
+%
+%   KERNEL.AXISSYMHELM2D('d', ZK) or KERNEL.AXISSYMHELM2D('double', ZK) 
+%   constructs the axissymmetric double-layer Helmholtz kernel with 
+%   wavenumber ZK.
+%
+%   KERNEL.AXISSYMHELM2D('sp', ZK) or KERNEL.AXISSYMHELM2D('sprime', ZK) 
+%   constructs the normal derivative of the axissymmetric single-layer 
+%   Helmholtz kernel with wavenumber ZK.
+%
+%   KERNEL.AXISSYMHELM2D('c', ZK, COEFS) or 
+%   KERNEL.AXISSYMHELM2D('combined', ZK, COEFS)
+%   constructs the combined-layer axissymmetric Helmholtz kernel with 
+%   wavenumber ZK and parameter COEFS, 
+%   i.e., COEFS(1)*KERNEL.AXISSYMHELM2D('d', ZK) + 
+%         COEFS(2)*KERNEL.AXISSYMHELM2D('s', ZK).
+%
+% NOTES: The axissymetric kernels are currently supported only for purely
+% real or purely imaginary wave numbers
+%
+% See also CHNK.AXISSYMHELM2D.KERN.
+
+if ( nargin < 1 )
+    error('Missing Laplace kernel type.');
+end
+
+obj = kernel();
+obj.name = 'axissymlaplace';
+obj.opdims = [1 1];
+
+switch lower(type)
+
+    case {'s', 'single'}
+        obj.type = 's';
+        obj.eval = @(s,t) chnk.axissymlap2d.kern(s, t, [0,0], 's', n);
+        obj.shifted_eval = @(s,t,o) chnk.axissymlap2d.kern(s, t, o, 's', n);
+        obj.fmm = [];
+        obj.sing = 'log';
+
+    case {'d', 'double'}
+        obj.type = 'd';
+        obj.eval = @(s,t) chnk.axissymlap2d.kern(s, t, [0,0], 'd', n);
+        obj.shifted_eval = @(s,t,o) chnk.axissymlap2d.kern(s, t, o, 'd', n);
+        obj.fmm = [];
+        obj.sing = 'log';
+
+    case {'sp', 'sprime'}
+        obj.type = 'sp';
+        obj.eval = @(s,t) chnk.axissymlap2d.kern(s, t, [0,0], 'sprime', n);
+        obj.shifted_eval = @(s,t,o) chnk.axissymlap2d.kern(s, t, o, 'sprime', n);
+        obj.fmm = [];
+        obj.sing = 'log';
+
+    case {'dp', 'dprime'}
+        obj.type = 'dp';
+        obj.eval = @(s,t) chnk.axissymlap2d.kern(s, t, [0,0], 'dprime', n);
+        obj.shifted_eval = @(s,t,o) chnk.axissymlap2d.kern(s, t, o, 'dprime', n);
+        obj.fmm = [];
+        obj.sing = 'hs';
+
+    case {'q', 'quad', 'quadruple', 'quadrupole'}
+        obj.type = 'q';
+        obj.eval = @(s,t) chnk.axissymlap2d.kern(s, t, [0,0], 'q', n);
+        obj.shifted_eval = @(s,t,o) chnk.axissymlap2d.kern(s, t, o, 'q', n);
+        obj.fmm = [];
+        obj.sing = 'hs';
+
+    case {'spp', 'sprimeprime'}
+        obj.type = 'spp';
+        obj.eval = @(s,t) chnk.axissymlap2d.kern(s, t, [0,0], 'spp', n);
+        obj.shifted_eval = @(s,t) chnk.axissymlap2d.kern(s, t, o, 'spp', n);
+        obj.fmm = [];
+        obj.sing = 'hs';
+
+    case {'q_sum_dp'}
+        obj.type = 'q_sum_dp';
+        obj.eval = @(s,t) chnk.axissymlap2d.kern(s, t, [0,0], 'q_sum_dp', n);
+        obj.shifted_eval = @(s,t,o) chnk.axissymlap2d.kern(s, t, o, 'q_sum_dp', n);
+        obj.fmm = [];
+        obj.sing = 'log';
+
+    case {'spp_sum_dp'}
+        obj.type = 'spp_sum_dp';
+        obj.eval = @(s,t) chnk.axissymlap2d.kern(s, t, [0,0], 'spp_sum_dp', n);
+        obj.shifted_eval = @(s,t) chnk.axissymlap2d.kern(s, t, o, 'spp_sum_dp', n);
+        obj.fmm = [];
+        obj.sing = 'log';
+
+
+    otherwise
+        error('Unknown axissym Laplace kernel type ''%s''.', type);
+
+end
+
+end
diff --git a/chunkie/+chnk/+axissymlap2d/gaus_agm.m b/chunkie/+chnk/+axissymlap2d/gaus_agm.m
new file mode 100644
index 00000000..dca93eff
--- /dev/null
+++ b/chunkie/+chnk/+axissymlap2d/gaus_agm.m
@@ -0,0 +1,38 @@
+function [rk0,re0] = gaus_agm(x)
+
+eps  = 1E-15;
+a    = sqrt(2./(x+1));
+delt = 1./sqrt(1-a.*a);
+aa0  = delt + sqrt(delt.*delt-1);
+bb0   = 1./(delt+sqrt(delt.*delt-1));
+a0  = ones(size(delt));
+b0  = 1./delt;
+
+
+fact = ((a0+b0)/2).^2;
+
+for i=1:1000
+    a1 = (a0+b0)/2;
+    b1 = sqrt(a0.*b0);
+
+    aa1 = (aa0+bb0)/2;
+    bb1 = sqrt(aa0.*bb0);
+    a0 = a1;
+    b0 = b1;
+    aa0 = aa1;
+    bb0 = bb1;
+
+    c0 = (a1-b1)/2;
+    fact = fact-(c0.*c0)*2^(i);
+    drel = abs(a0-b0)./abs(a0);
+    drel2= abs(aa0-bb0)./abs(aa0);
+    if (max(drel+drel2) <2*eps)
+        break
+    end
+
+end
+
+rk0 = pi./(2*aa0.*sqrt(1-a.*a));
+re0 = pi*fact./(2*a0);
+
+end
\ No newline at end of file
diff --git a/chunkie/+chnk/+axissymlap2d/gfunc.m b/chunkie/+chnk/+axissymlap2d/gfunc.m
new file mode 100644
index 00000000..e89acabd
--- /dev/null
+++ b/chunkie/+chnk/+axissymlap2d/gfunc.m
@@ -0,0 +1,64 @@
+function [gval,gdz,gdr,gdrp,gdzz,gdrrp,gdrz,gdrpz] = gfunc(r,rp,dr,dz,n)
+    % chnk.axissymlap2d.gfunc
+    % green's function kernel normalized by area
+    domega = 1/(2*pi);
+    %wn = @(n) 2*pi^(n/2)/gamma(n/2);
+    t0 = (r == 0);
+    s0 = (rp == 0);
+    st0 = ~((~t0).*(~s0));
+    dz2 = dz.^2;
+    gval0 = 1/2*rp./(sqrt(r.^2+rp.^2+dz2));
+    gdrp0 = 1/2*rp.^2./(sqrt(rp.^2+dz.^2).^3); 
+    gdz0  = 1/2*rp.*dz./(sqrt(rp.^2+dz.^2).^3);
+
+    t = (dz.^2+dr.^2)./(2.*r.*rp);
+    factor = (rp./r).^((n-2)/2);
+    %[q0,~,q0d,q0dd] = chnk.axissymlap2d.qleg_half(t); % TODO: q0dd for D'
+    %[~,~,~,q0dd] = chnk.axissymlap2d.qleg_half(t); % TODO: q0dd for D'
+    inear = (t<1e-2);
+    xm1n = t(inear);
+    xm1f = t(~inear);
+    [q0f,~,qdf,qddf] = chnk.axissymlap2d.runbackward(xm1f,n); 
+    [q0n,~,qdn,qddn] = chnk.axissymlap2d.runforward(xm1n,n); 
+    q0 = zeros(size(t)); % q0 = Q_{n/2-2}
+    q0d = q0; q0dd = q0;
+    q0(inear) = q0n; q0(~inear) = q0f;
+    q0d(inear) = qdn; q0d(~inear) = qdf;
+    q0dd(inear) = qddn; q0dd(~inear) = qddf;
+
+    gval = domega*factor.*q0;
+    gdz  =-domega*factor.*q0d ...
+           ./(rp.*r).*dz;
+    rpfac = -r*(n-2)/2.*q0+(-(1+t).*r+rp).*q0d;
+    gdrp  = -domega*factor./(rp.*r) ...
+          .*rpfac;
+    rfac = -rp*(n-2)/2.*q0+(-(1+t).*rp+r).*q0d;
+    gdr  = -domega*factor./(rp.*r) ...
+          .*rfac;
+
+    gval(st0) = gval0(st0);
+    gdr(st0) = 0;
+    gdrp(st0) = gdrp0(st0);
+    gdz(st0) = gdz0(st0);
+     
+    % t derivatives * rrp
+    tdr = r-rp.*(1+t);
+    tdrp = rp-(1+t).*r;
+    tdz = -dz;
+    rrp = rp.*r;
+
+    % compute for hessian
+    gdzz = domega*factor.*(q0dd.*(dz./rrp).^2+q0d./rrp); 
+
+    %q0dterm = 2*(1+t)-(r./rp+rp./r)*3/2;
+    q0dterm = (n-1)*(1+t)-(r./rp+rp./r)*n/2;
+    gdrrp = (n-2)^2*q0/4 + q0d.*q0dterm + q0dd./rrp.*tdr.*tdrp;
+    gdrrp = -gdrrp*domega.*factor./rrp;
+    
+    gdrz = q0dd.*tdr - q0d.*rp*n/2;
+    gdrz = -gdrz*domega.*factor.*dz./rrp.^2; 
+
+    gdrpz = q0dd.*tdrp - q0d.*r*n/2;
+    gdrpz = -gdrpz*domega.*factor.*dz./rrp.^2; 
+    
+end
\ No newline at end of file
diff --git a/chunkie/+chnk/+axissymlap2d/green.m b/chunkie/+chnk/+axissymlap2d/green.m
new file mode 100644
index 00000000..326ae03c
--- /dev/null
+++ b/chunkie/+chnk/+axissymlap2d/green.m
@@ -0,0 +1,57 @@
+function [val, grad, hess] = green(src, targ, origin, n)
+%CHNK.AXISSYMLAP2D.GREEN evaluate the Laplace green's function
+% for the given sources and targets. 
+%
+% Note: that the first coordinate is r, and the second z.
+% The code relies on precomputed tables and hence loops are required for 
+% computing various pairwise interactions.
+% Finally, the code is not efficient in the sense that val, grad, hess 
+% are always internally computed independent of nargout
+%
+% Since the kernels are not translationally invariant in r, the size
+% of the gradient is 3, for d/dr, d/dr', and d/dz
+%
+% Similarly the hessian is of size 6 and ordered as 
+%   d_{rr}, d_{r'r'}, d_{zz}, d_{rr'}, d_{rz}, d_{r'z}
+
+[~, ns] = size(src);
+[~, nt] = size(targ);
+
+vtmp = zeros(nt, ns);
+gtmp = zeros(nt, ns, 3);
+htmp = zeros(nt, ns, 6);
+
+rt = repmat(targ(1,:).',1,ns);
+rs = repmat(src(1,:),nt,1);
+dz = repmat(src(2,:),nt,1)-repmat(targ(2,:).',1,ns);
+r  = (rt + origin(1));
+rp = (rs + origin(1));
+dr = (rs-rt);
+
+[gval,gdz,gdr,gdrp,gdzz,gdrrp,gdrz,gdrpz] = chnk.axissymlap2d.gfunc(r,rp,dr,dz,n);
+
+vtmp = gval;
+gtmp(:,:,1) = gdr;
+gtmp(:,:,2) = gdrp;
+gtmp(:,:,3) = gdz;
+
+htmp(:,:,3) = gdzz;
+htmp(:,:,4) = gdrrp;
+htmp(:,:,5) = gdrz;
+htmp(:,:,6) = gdrpz;
+
+if nargout > 0
+    val = vtmp;
+end
+
+if nargout > 1
+    grad = gtmp;
+end
+
+if nargout > 2
+    hess = htmp;
+end
+
+end
+
+
diff --git a/chunkie/+chnk/+axissymlap2d/kern.m b/chunkie/+chnk/+axissymlap2d/kern.m
new file mode 100644
index 00000000..617b1a42
--- /dev/null
+++ b/chunkie/+chnk/+axissymlap2d/kern.m
@@ -0,0 +1,137 @@
+function submat = kern(srcinfo, targinfo, origin, type, dimension)
+%CHNK.AXISSYMLAP2D.KERN axissymmetric Laplace layer potential kernels in 2D
+% 
+% Syntax: submat = chnk.axissymlap2d.kern(srcinfo,targingo,type)
+%
+% Let x be targets and y be sources for these formulas, with
+% n_x and n_y the corresponding unit normals at those points
+% (if defined). Note that the normal information is obtained
+% by taking the perpendicular to the provided tangential deriviative
+% info and normalizing  
+%
+% Here the first and second components correspond to the r and z
+% coordinates respectively. 
+%
+% Kernels based on G(x,y) = \int_{0}^{\pi} 1/(d(t)) \, dt \, 
+% where d(t) = \sqrt(r^2 + r'^2 - 2rr' \cos(t) + (z-z')^2) with
+% x = (r,z), and y = (r',z')
+%
+% D(x,y) = \nabla_{n_y} G(x,y)
+% S(x,y) = G(x,y)
+% S'(x,y) = \nabla_{n_x} G(x,y)
+% D'(x,y) = \nabla_{n_x} \nabla_{n_y} G(x,y)
+%
+% Input:
+%   srcinfo - description of sources in ptinfo struct format, i.e.
+%                ptinfo.r - positions (2,:) array
+%                ptinfo.d - first derivative in underlying
+%                     parameterization (2,:)
+%                ptinfo.d2 - second derivative in underlying
+%                     parameterization (2,:)
+%   targinfo - description of targets in ptinfo struct format,
+%                if info not relevant (d/d2) it doesn't need to
+%                be provided. sprime requires tangent info in
+%                targinfo.d
+%   type - string, determines kernel type
+%                type == 'd', double layer kernel D
+%                type == 's', single layer kernel S
+%                type == 'sprime', normal derivative of single
+%                      layer S'
+% Output:
+%   submat - the evaluation of the selected kernel for the
+%            provided sources and targets. the number of
+%            rows equals the number of targets and the
+%            number of columns equals the number of sources  
+src = srcinfo.r; 
+targ = targinfo.r;
+[~, ns] = size(src);
+[~, nt] = size(targ);
+
+if strcmpi(type, 'd')
+    srcnorm = srcinfo.n;
+    [~, grad] = chnk.axissymlap2d.green(src, targ, origin, dimension);
+    nx = repmat(srcnorm(1,:), nt, 1);
+    ny = repmat(srcnorm(2,:), nt, 1);
+    % Due to lack of translation invariance in r, no sign flip needed, 
+    % as gradient is computed with repsect to r'
+    submat = (grad(:,:,2).*nx + grad(:,:,3).*ny);
+end
+
+if strcmpi(type, 'dprime')
+    targnorm = targinfo.n;
+    srcnorm = srcinfo.n;
+    [~,~,hess] = chnk.axissymlap2d.green(src, targ, origin, dimension);
+    nxsrc = repmat(srcnorm(1,:),nt,1);
+    nysrc = repmat(srcnorm(2,:),nt,1);
+    nxtarg = repmat((targnorm(1,:)).',1,ns);
+    nytarg = repmat((targnorm(2,:)).',1,ns);
+    submat = hess(:,:,4).*nxsrc.*nxtarg + hess(:,:,5).*nysrc.*nxtarg ... 
+              - hess(:,:,6).*nxsrc.*nytarg + hess(:,:,3).*nysrc.*nytarg;
+end
+
+if strcmpi(type, 'sprime')
+    targnorm = targinfo.n;
+    [~, grad] = chnk.axissymlap2d.green(src, targ, origin, dimension);
+    nx = repmat((targnorm(1,:)).',1,ns);
+    ny = repmat((targnorm(2,:)).',1,ns);
+    submat = (grad(:,:,1).*nx - grad(:,:,3).*ny);
+end
+
+if strcmpi(type, 's')
+    submat = chnk.axissymlap2d.green(src, targ, origin, dimension);  
+end
+
+if strcmpi(type,'sprimeprime')
+    targnorm = targinfo.n(:,:);
+    [~,~,hess] = chnk.axissymlap2d.green(src, targ, origin, dimension);
+    nxtarg = repmat((targnorm(1,:)).',1,ns);
+    nytarg = repmat((targnorm(2,:)).',1,ns);
+    submat = hess(:,:,1).*nxtarg.*nxtarg - 2*hess(:,:,5).*nytarg.*nxtarg ...
+                                          + hess(:,:,3).*nytarg.*nytarg;
+end
+
+if strcmpi(type,'q')
+    srcnorm = srcinfo.n(:,:);
+    [~,~,hess] = chnk.axissymlap2d.green(src, targ, origin, dimension);
+    nxsrc = repmat(srcnorm(1,:),nt,1);
+    nysrc = repmat(srcnorm(2,:),nt,1);
+    submat = hess(:,:,2).*nxsrc.*nxsrc + 2*hess(:,:,6).*nysrc.*nxsrc ... 
+                                      + hess(:,:,3).*nysrc.*nysrc;
+end
+
+if strcmpi(type,'q_sum_dp')
+    % D'
+    targnorm = targinfo.n;
+    srcnorm = srcinfo.n;
+    [~,~,hess] = chnk.axissymlap2d.green(src, targ, origin, dimension);
+    nxsrc = repmat(srcnorm(1,:),nt,1);
+    nysrc = repmat(srcnorm(2,:),nt,1);
+    nxtarg = repmat((targnorm(1,:)).',1,ns);
+    nytarg = repmat((targnorm(2,:)).',1,ns);
+    submat = hess(:,:,4).*nxsrc.*nxtarg + hess(:,:,5).*nysrc.*nxtarg ... 
+              - hess(:,:,6).*nxsrc.*nytarg + hess(:,:,3).*nysrc.*nytarg;
+
+    % add Q
+    submat = submat + (hess(:,:,2).*nxsrc.*nxsrc + 2*hess(:,:,6).*nysrc.*nxsrc ... 
+                                      + hess(:,:,3).*nysrc.*nysrc);
+end
+
+if strcmpi(type,'spp_sum_dp')
+    % D'
+    targnorm = targinfo.n;
+    srcnorm = srcinfo.n;
+    [~,~,hess] = chnk.axissymlap2d.green(src, targ, origin, dimension);
+    nxsrc = repmat(srcnorm(1,:),nt,1);
+    nysrc = repmat(srcnorm(2,:),nt,1);
+    nxtarg = repmat((targnorm(1,:)).',1,ns);
+    nytarg = repmat((targnorm(2,:)).',1,ns);
+    submat = hess(:,:,4).*nxsrc.*nxtarg + hess(:,:,5).*nysrc.*nxtarg ... 
+              - hess(:,:,6).*nxsrc.*nytarg + hess(:,:,3).*nysrc.*nytarg;
+
+    % add S''
+    submat = submat + (hess(:,:,1).*nxtarg.*nxtarg + 2*hess(:,:,5).*nytarg.*nxtarg ...
+                                          + hess(:,:,3).*nytarg.*nytarg);
+end
+
+
+end
\ No newline at end of file
diff --git a/chunkie/+chnk/+axissymlap2d/qeval0_far.m b/chunkie/+chnk/+axissymlap2d/qeval0_far.m
new file mode 100644
index 00000000..3c3727cd
--- /dev/null
+++ b/chunkie/+chnk/+axissymlap2d/qeval0_far.m
@@ -0,0 +1,24 @@
+function [q0] = qeval0_far(x)
+
+        coefs = ...
+       [2.2214414690791831235E0,
+        0,                    
+        0.41652027545234683566E0,                    
+        0,                    
+        0.22778452563800217575E0,                    
+        0,                    
+        0.15660186137612649583E0,                    
+        0,                    
+        0.11928657409509635424E0]; 
+
+
+        t = 1./x;
+        t0 = 1./sqrt(x);
+        q0 = zeros(size(t));
+
+        for i=1:9
+                q0 = q0 + t0*coefs(i);
+                t0 = t0.*t;        
+        end
+
+end
\ No newline at end of file
diff --git a/chunkie/+chnk/+axissymlap2d/qeval0_near.m b/chunkie/+chnk/+axissymlap2d/qeval0_near.m
new file mode 100644
index 00000000..185d47e8
--- /dev/null
+++ b/chunkie/+chnk/+axissymlap2d/qeval0_near.m
@@ -0,0 +1,34 @@
+function [q0] = qeval0_near(t)
+
+        coefs = ...
+             [1.7328679513998632735E0,
+            -0.5000000000000000000E0,
+           -0.09160849392498290919E0,
+           0.062500000000000000000E0,
+           0.019905513916401443210E0,
+          -0.017578125000000000000E0,
+          -0.006097834693194945559E0,
+           0.006103515625000000000E0,
+          0.0021674343381175963469E0,
+         -0.0023365020751953125000E0,
+         -0.0008357538695841108955E0,
+          0.0009462833404541015625E0,
+          0.0003390851133264318725E0,
+        -0.00039757043123245239258E0,
+        -0.00014242013770157621830E0,
+         0.00017140153795480728149E0,
+         0.00006133159221782055041E0,
+        -0.00007532294148404616863E0];
+
+
+        b = log(t);
+        v = 1;
+        
+        q0 = zeros(size(t));
+
+        for i=1:9
+                q0 = q0 + v*coefs(2*i-1)+v.*b*coefs(2*i);
+                v  = v.*t;        
+        end
+
+end
\ No newline at end of file
diff --git a/chunkie/+chnk/+axissymlap2d/qeval0der_near.m b/chunkie/+chnk/+axissymlap2d/qeval0der_near.m
new file mode 100644
index 00000000..0ab62544
--- /dev/null
+++ b/chunkie/+chnk/+axissymlap2d/qeval0der_near.m
@@ -0,0 +1,36 @@
+function [q0d] = qeval0der_near(t)
+
+        coefs = ...
+             [1.7328679513998632735d0,
+            -0.5000000000000000000d0,
+           -0.09160849392498290919d0,
+           0.062500000000000000000d0,
+           0.019905513916401443210d0,
+          -0.017578125000000000000d0,
+          -0.006097834693194945559d0,
+           0.006103515625000000000d0,
+          0.0021674343381175963469d0,
+         -0.0023365020751953125000d0,
+         -0.0008357538695841108955d0,
+          0.0009462833404541015625d0,
+          0.0003390851133264318725d0,
+        -0.00039757043123245239258d0,
+        -0.00014242013770157621830d0,
+         0.00017140153795480728149d0,
+         0.00006133159221782055041d0,
+        -0.00007532294148404616863d0];
+
+        b = log(t);
+        v = 1;
+
+        q0d = coefs(2)./t;
+
+        for i=2:9
+                q0d = q0d + (i-1)* ...
+                   (v*coefs(2*i-1)+v.*b*coefs(2*i)) ...
+                       +v*coefs(2*i);
+                v  = v.*t;        
+        end
+     
+
+end
\ No newline at end of file
diff --git a/chunkie/+chnk/+axissymlap2d/qeval1_far.m b/chunkie/+chnk/+axissymlap2d/qeval1_far.m
new file mode 100644
index 00000000..27286bcc
--- /dev/null
+++ b/chunkie/+chnk/+axissymlap2d/qeval1_far.m
@@ -0,0 +1,24 @@
+function [q1] = qeval1_far(x)
+
+        coefs = ...
+       [0.55536036726979578088E0,                    
+        0,                    
+        0.26032517215771677229E0,                    
+        0,                    
+        0.17083839422850163181E0,                    
+        0,                    
+        0.12723901236810277786E0,                    
+        0,                    
+        0.10139358798083190111E0];  
+
+
+        t = 1./x;
+        t0 = 1./sqrt(x).^3;
+        q1 = zeros(size(t));
+
+        for i=1:9
+                q1 = q1 + t0*coefs(i);
+                t0 = t0.*t;        
+        end
+
+end
\ No newline at end of file
diff --git a/chunkie/+chnk/+axissymlap2d/qeval1_near.m b/chunkie/+chnk/+axissymlap2d/qeval1_near.m
new file mode 100644
index 00000000..cd506473
--- /dev/null
+++ b/chunkie/+chnk/+axissymlap2d/qeval1_near.m
@@ -0,0 +1,34 @@
+function [q1] = qeval1_near(t)
+
+        coefs = ...
+            [-0.2671320486001367265d0,
+            -0.5000000000000000000d0,
+             0.5248254817749487276d0,
+           -0.18750000000000000000d0,
+           -0.04098835652733573868d0,
+           0.029296875000000000000d0, 
+           0.009513531070472923783d0,
+          -0.008544921875000000000d0,
+         -0.0029774361551467310174d0, 
+          0.0030040740966796875000d0, 
+          0.0010682069932178195667d0,
+         -0.0011565685272216796875d0,
+         -0.0004138797762860294822d0, 
+         0.00046985596418380737305d0, 
+          0.00016838776925222835322d0,
+        -0.00019777100533246994019d0,
+        -0.00007084821236213522235d0, 
+         0.00008536600034858565778d0];
+
+
+        b = log(t);
+        v = 1;
+        
+        q1 = zeros(size(t));
+
+        for i=1:9
+                q1 = q1 + v*coefs(2*i-1)+v.*b*coefs(2*i);
+                v  = v.*t;        
+        end
+
+end
\ No newline at end of file
diff --git a/chunkie/+chnk/+axissymlap2d/qget_zero.m b/chunkie/+chnk/+axissymlap2d/qget_zero.m
new file mode 100644
index 00000000..0e56a581
--- /dev/null
+++ b/chunkie/+chnk/+axissymlap2d/qget_zero.m
@@ -0,0 +1,16 @@
+function [qm,qmp] = qget_zero(x)
+
+        a = sqrt(2./(x+1));
+        [fF,fE] = chnk.axissymlap2d.gaus_agm(x);
+
+        q0 = a.*fF;
+        q1 = x.*q0-2*fE./(a);
+
+        qa = q0;
+        qb = q1;
+
+        qmm = 0;
+        qm  = q0;
+        qmp = q1;
+
+end
\ No newline at end of file
diff --git a/chunkie/+chnk/+axissymlap2d/qleg_half.m b/chunkie/+chnk/+axissymlap2d/qleg_half.m
new file mode 100644
index 00000000..d1fa2bd4
--- /dev/null
+++ b/chunkie/+chnk/+axissymlap2d/qleg_half.m
@@ -0,0 +1,39 @@
+function [q0,q1,q0d,q0dd] = qleg_half(t)
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%
+%   caveat utilitor: this function evaluates Q_{-1/2} and 
+%                    Q_{1/2} at t+1;
+%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+    inear = find(t<0.01);
+    ifar  = find(t>100);
+    imid  = find((t>=0.01).*(t<=100));
+    [qmm,qmpm] = chnk.axissymlap2d.qget_zero(1+t(imid));
+    [qmn] = chnk.axissymlap2d.qeval0_near(t(inear));
+    [qmpn] = chnk.axissymlap2d.qeval1_near(t(inear));
+    [qmdn] = chnk.axissymlap2d.qeval0der_near(t(inear));
+    [qmf] = chnk.axissymlap2d.qeval0_far(1+t(ifar));
+    [qmpf] = chnk.axissymlap2d.qeval1_far(1+t(ifar));
+
+
+    
+    q0 = zeros(size(t));
+    q0(imid)  = qmm;
+    q0(inear) = qmn;
+    q0(ifar)  = qmf;
+
+    q1 = zeros(size(t));
+    q1(imid)  = qmpm;
+    q1(inear) = qmpn;
+    q1(ifar)  = qmpf;
+
+    q0d = -(1/2)*q1 + (1/2)*(1+t).*q0;
+    q0d = q0d./(-t.*(t+2));
+    q0d(inear) = qmdn;
+
+    q1d = (1/2)*q0 - (1/2)*(1+t).*q1;
+    q1d = q1d./(-t.*(t+2));
+    q0dd = -(1/2)*q1d + (1/2)*q0 + (5/2)*(1+t).*q0d;
+    q0dd = q0dd./(-t.*(t+2));
+end
\ No newline at end of file
diff --git a/chunkie/+chnk/+axissymlap2d/runbackward.m b/chunkie/+chnk/+axissymlap2d/runbackward.m
new file mode 100644
index 00000000..21388767
--- /dev/null
+++ b/chunkie/+chnk/+axissymlap2d/runbackward.m
@@ -0,0 +1,70 @@
+function [q0,q1,qd0,qdd0] = runbackward(xm1,n)
+    % returns q0 = Q_{n/2-2}
+    if rem(n, 2) == 0
+        % n is even
+        first_term = 1/2*log((xm1+2)./xm1); % Q0
+        m = n/2-1;
+    else
+        % n is odd, n>1
+        first_term = chnk.axissymlap2d.qleg_half(xm1); % Q_{-1/2}
+        m = (n-1)/2;
+    end
+    j = 0;
+    jmax = 300;
+    tol = 1e-8; 
+    t = repmat(zeros(size(xm1)),1,1,2);
+    t(:,:,1) = ones(size(xm1));
+    t(:,:,2) = zeros(size(xm1));
+    onezero = t;
+    told = repmat(zeros(size(xm1)),1,1,2);
+    told(:,:,1) = realmax*ones(size(xm1));
+    told(:,:,2) = ones(size(xm1));
+    x = xm1+1;
+    while j<jmax
+        v = n/2-1+j; %m+j-1/2;
+        qv1 = t(:,:,2);
+        qv0 = t(:,:,1); 
+        qv2 = (2*v+3)/(v+2)*x.*qv1-(v+1)/(v+2)*qv0;
+        t(:,:,1) = qv1;
+        t(:,:,2) = qv2;
+        vecn = vecnorm(t,2,3);
+        t(:,:,1) = t(:,:,1)./vecn;
+        t(:,:,2) = t(:,:,2)./vecn;
+        t = t./t(:,:,2); 
+        if max(abs(told(:,:,1)-t(:,:,1)),[],'all')<tol
+            break;
+        end 
+        told = t;
+        j = j+1;
+    end
+    t = onezero;
+    for ii = m+j:-1:m
+        qv1 = t(:,:,1);
+        qv2 = t(:,:,2);
+        qv0 = (2*v+3)/(v+1)*x.*qv1-(v+2)/(v+1)*qv2;
+        t(:,:,1) = qv0;
+        t(:,:,2) = qv1;
+        t = t./qv0;
+        v = v-1;
+    end
+    qm1 = t(:,:,2);
+    for ii = m-1:-1:1
+        qv1 = t(:,:,1);
+        qv2 = t(:,:,2);
+        qv0 = (2*v+3)/(v+1)*x.*qv1-(v+2)/(v+1)*qv2;
+        t(:,:,1) = qv0;
+        t(:,:,2) = qv1;
+        v = v-1;
+    end
+    
+    scale = first_term./t(:,:,1);
+    q0 = scale; 
+    q1 = qm1.*scale; 
+    v = n/2-2; %
+    qd0 = (-v-1)*q1+(v+1)*(1+xm1).*q0; 
+    qd0 = -qd0./(2+xm1)./xm1;
+    qd1 = (v+1)*q0-(v+1)*(1+xm1).*q1;
+    qd1 = -qd1./(2+xm1)./xm1;
+    qdd0 = (-v-1)*qd1+(v+1)*q0+(v+3)*(1+xm1).*qd0;
+    qdd0 = -qdd0./(2+xm1)./xm1;
+end
\ No newline at end of file
diff --git a/chunkie/+chnk/+axissymlap2d/runforward.m b/chunkie/+chnk/+axissymlap2d/runforward.m
new file mode 100644
index 00000000..03b8af02
--- /dev/null
+++ b/chunkie/+chnk/+axissymlap2d/runforward.m
@@ -0,0 +1,33 @@
+function [q0,q1,qd0,qdd0] = runforward(xm1,n)
+    % returns q0 = Q_{n/2-2}
+    if rem(n, 2) == 0
+        % n is even
+        Q1 = 1/2*log((xm1+2)./xm1); % Q0
+        Q2 = (xm1+1)/2.*log((xm1+2)./xm1)-1; % Q1
+        m = n/2-1;
+        v = 0;
+    else
+        % n is odd, n>1
+        [Q1,Q2,~,~] = chnk.axissymlap2d.qleg_half(xm1); % Q_{-1/2}
+        m = (n-1)/2;
+        v = -1/2;
+    end
+    t(:,:,1) = Q1;
+    t(:,:,2) = Q2;
+    while v < m-3/2
+        qv0 = t(:,:,1);
+        qv1 = t(:,:,2);
+        qv2 = (2*v+3)/(v+2)*(xm1+1).*qv1-(v+1)/(v+2)*qv0;
+        t(:,:,1) = qv1; 
+        t(:,:,2) = qv2;
+        v = v+1;
+    end
+    q0 = t(:,:,1); % Q_{m-3/2}
+    q1 = t(:,:,2); % Q_{m-1/2}
+    qd0 = (-v-1)*q1+(v+1)*(1+xm1).*q0; 
+    qd0 = -qd0./(2+xm1)./xm1;
+    qd1 = (v+1)*q0-(v+1)*(1+xm1).*q1;
+    qd1 = -qd1./(2+xm1)./xm1;
+    qdd0 = (-v-1)*qd1+(v+1)*q0+(v+3)*(1+xm1).*qd0;
+    qdd0 = -qdd0./(2+xm1)./xm1;
+end
diff --git a/chunkie/@kernel/axissymlap2d.m b/chunkie/@kernel/axissymlap2d.m
new file mode 100644
index 00000000..af1d6487
--- /dev/null
+++ b/chunkie/@kernel/axissymlap2d.m
@@ -0,0 +1,98 @@
+function obj = axissymlap2d(type, n)
+%KERNEL.AXISSYMHELM2D   Construct the axissymmetric Helmholtz kernel.
+%   KERNEL.AXISSYMHELM2D('s', ZK) or KERNEL.AXISSYMHELM2D('single', ZK) 
+%   constructs the axissymmetric single-layer Helmholtz kernel with 
+%   wavenumber ZK.
+%
+%   KERNEL.AXISSYMHELM2D('d', ZK) or KERNEL.AXISSYMHELM2D('double', ZK) 
+%   constructs the axissymmetric double-layer Helmholtz kernel with 
+%   wavenumber ZK.
+%
+%   KERNEL.AXISSYMHELM2D('sp', ZK) or KERNEL.AXISSYMHELM2D('sprime', ZK) 
+%   constructs the normal derivative of the axissymmetric single-layer 
+%   Helmholtz kernel with wavenumber ZK.
+%
+%   KERNEL.AXISSYMHELM2D('c', ZK, COEFS) or 
+%   KERNEL.AXISSYMHELM2D('combined', ZK, COEFS)
+%   constructs the combined-layer axissymmetric Helmholtz kernel with 
+%   wavenumber ZK and parameter COEFS, 
+%   i.e., COEFS(1)*KERNEL.AXISSYMHELM2D('d', ZK) + 
+%         COEFS(2)*KERNEL.AXISSYMHELM2D('s', ZK).
+%
+% NOTES: The axissymetric kernels are currently supported only for purely
+% real or purely imaginary wave numbers
+%
+% See also CHNK.AXISSYMHELM2D.KERN.
+
+if ( nargin < 1 )
+    error('Missing Laplace kernel type.');
+end
+
+obj = kernel();
+obj.name = 'axissymlaplace';
+obj.opdims = [1 1];
+
+switch lower(type)
+
+    case {'s', 'single'}
+        obj.type = 's';
+        obj.eval = @(s,t) chnk.axissymlap2d.kern(s, t, [0,0], 's', n);
+        obj.shifted_eval = @(s,t,o) chnk.axissymlap2d.kern(s, t, o, 's', n);
+        obj.fmm = [];
+        obj.sing = 'log';
+
+    case {'d', 'double'}
+        obj.type = 'd';
+        obj.eval = @(s,t) chnk.axissymlap2d.kern(s, t, [0,0], 'd', n);
+        obj.shifted_eval = @(s,t,o) chnk.axissymlap2d.kern(s, t, o, 'd', n);
+        obj.fmm = [];
+        obj.sing = 'log';
+
+    case {'sp', 'sprime'}
+        obj.type = 'sp';
+        obj.eval = @(s,t) chnk.axissymlap2d.kern(s, t, [0,0], 'sprime', n);
+        obj.shifted_eval = @(s,t,o) chnk.axissymlap2d.kern(s, t, o, 'sprime', n);
+        obj.fmm = [];
+        obj.sing = 'log';
+
+    case {'dp', 'dprime'}
+        obj.type = 'dp';
+        obj.eval = @(s,t) chnk.axissymlap2d.kern(s, t, [0,0], 'dprime', n);
+        obj.shifted_eval = @(s,t,o) chnk.axissymlap2d.kern(s, t, o, 'dprime', n);
+        obj.fmm = [];
+        obj.sing = 'hs';
+
+    case {'spp', 'sprimeprime'}
+        obj.type = 'spp';
+        obj.eval = @(s,t) chnk.axissymlap2d.kern(s, t, [0,0], 'sprimeprime', n);
+        obj.shifted_eval = @(s,t,o) chnk.axissymlap2d.kern(s, t, o, 'sprimeprime', n);
+        obj.fmm = [];
+        obj.sing = 'hs';
+
+    case {'q'}
+        obj.type = 'q';
+        obj.eval = @(s,t) chnk.axissymlap2d.kern(s, t, [0,0], 'q', n);
+        obj.shifted_eval = @(s,t,o) chnk.axissymlap2d.kern(s, t, o, 'q', n);
+        obj.fmm = [];
+        obj.sing = 'hs';
+
+    case {'q_sum_dp'}
+        obj.type = 'q_sum_dp';
+        obj.eval = @(s,t) chnk.axissymlap2d.kern(s, t, [0,0], 'q_sum_dp', n);
+        obj.shifted_eval = @(s,t,o) chnk.axissymlap2d.kern(s, t, o, 'q_sum_dp', n);
+        obj.fmm = [];
+        obj.sing = 'hs';
+
+    case {'spp_sum_dp'}
+        obj.type = 'spp_sum_dp';
+        obj.eval = @(s,t) chnk.axissymlap2d.kern(s, t, [0,0], 'spp_sum_dp', n);
+        obj.shifted_eval = @(s,t,o) chnk.axissymlap2d.kern(s, t, o, 'spp_sum_dp', n);
+        obj.fmm = [];
+        obj.sing = 'hs';
+     
+    otherwise
+        error('Unknown axissym Laplace kernel type ''%s''.', type);
+
+end
+
+end
diff --git a/chunkie/@kernel/kernel.m b/chunkie/@kernel/kernel.m
index 39efb846..911644cc 100644
--- a/chunkie/@kernel/kernel.m
+++ b/chunkie/@kernel/kernel.m
@@ -113,8 +113,10 @@
                       obj = kernel.zeros(varargin{:});
                   case {'nans', 'nan'}
                       obj = kernel.nans(varargin{:});
-                  case {'axis sym helmholtz', 'axissymh', 'axissymhelm'}
+                  case {'axis sym helmholtz', 'axissymh', 'axissymhelm', 'axissymhelmholtz'}
                       obj = kernel.axissymhelm2d(varargin{:});
+                  case {'axis sym laplace', 'axissyml', 'axissymlap', 'axissymlaplace'}
+                      obj = kernel.axissymlap2d(varargin{:});
                   case {'axis sym helmholtz difference', 'axissymhdiff' ...
                            'axissymhelmdiff', 'axissymhelm_diff'}
                       obj = kernel.axissymhelm2ddiff(varargin{:});   
@@ -165,6 +167,7 @@
         obj = stok2d(varargin);
         obj = elast2d(varargin);
         obj = axissymhelm2d(varargin);
+        obj = axissymlap2d(varargin);
         obj = axissymhelm2ddiff(varargin);
         obj = helm2dquas(varargin);
         obj = zeros(varargin);
diff --git a/chunkie/demo/demo_axissymlap_analytic.m b/chunkie/demo/demo_axissymlap_analytic.m
new file mode 100644
index 00000000..a503d033
--- /dev/null
+++ b/chunkie/demo/demo_axissymlap_analytic.m
@@ -0,0 +1,121 @@
+%clearvars; clc;
+iftorus = 0;
+cparams = [];
+cparams.eps = 1.0e-10;
+cparams.nover = 1;
+if ~iftorus % sphere
+    cparams.ta = -pi/2;
+    cparams.tb = pi/2;
+    center = [0;2];
+    cparams.ifclosed = false;
+else % torus
+    cparams.ta = 0;
+    cparams.tb = 2*pi;
+    center = [3;0];
+    cparams.ifclosed = true;
+end
+cparams.maxchunklen = 0.5;
+radius = 1;
+fcurve = @(t) radius*[cos(t(:).'); sin(t(:).')];
+
+chnkrhalf = chunkerfunc(fcurve,cparams);
+chnkrhalf = chnkrhalf.move(-center); 
+
+figure(1);clf;plot(chnkrhalf);hold on;plot(chnkrhalf,'bo');quiver(chnkrhalf,'y');
+
+% analytic solution in n-dimension
+ndim = 5;
+ftrue = @(s) s.r(1,:).^2-(ndim-1)*s.r(2,:).^2; % analytical solution
+fdn = @(s) -2*(ndim-1)*s.r(2,:); % (2*x1,2*x2,...,-2*(ndim-1)*z)*(0,0,1) = -2*(ndim-1)*z
+fdnn = @(s) -2*(ndim-1)*ones(size(s.r(2,:)));
+
+rhs = 2*ftrue(chnkrhalf).';
+%kerns = kernel('axissymlaplace','s',ndim);
+kerndp = kernel('axissymlaplace','dp',ndim);
+kernd = kernel('axissymlaplace','d',ndim);
+D1 = 2*kernd;
+mat = chunkermat(chnkrhalf,D1);
+A = mat+eye(chnkrhalf.npt);
+sig = gmres(A, rhs, [], 1e-14, 200);
+% place targets
+x1 = linspace(center(1),center(1)+radius,200);
+x2 = linspace(center(2)-radius,center(2)+radius,200);
+[xx,yy] = meshgrid(x1,x2);
+t = [];
+t.r = [xx(:).'; yy(:).'];
+in = chunkerinterior(chnkrhalf,{x1,x2});
+t.r = t.r(:,in(:));
+ucomp = kernd.eval(chnkrhalf,t)*(sig.*chnkrhalf.wts(:));
+utrue = ftrue(t).';
+err = log10(abs((ucomp-utrue)./max(abs(utrue))));
+
+%% calculate derivative via finite difference
+eps = 1e-8;
+te = [];
+te.r = [xx(:).'; yy(:).']+[0;1]*eps;
+te.r = te.r(:,in(:));
+te.n = zeros(size(t.r))+[0;1];
+ucompe = kerndp.eval(chnkrhalf,te)*(sig.*chnkrhalf.wts(:));
+dudn_fd = (ucompe-ucomp)./eps;
+%% true derivative
+dudntrue = fdn(t).';
+errdudn_fd = log10(abs((dudntrue-dudn_fd)./dudntrue));
+
+%% calculate derivative via Dp inside domain
+t.n = zeros(size(t.r))+[0;1];
+upcompin = kerndp.eval(chnkrhalf,t)*(sig.*chnkrhalf.wts(:));
+errdudn = log10(abs((dudntrue-upcompin)./dudntrue));
+
+%% calculate derivative via Dp on boundary
+%uptrue = [2,-4]*(chnkrhalf.r(:,:).*chnkrhalf.n(:,:));
+%targpt = []; targpt.r = [0;1]; targpt.n = [0;1];
+%opts = []; opts.sing = 'hs';
+%mat = chunkermat(chnkrhalf,Dp,opts);
+%upcomp = mat*(sig.*chnkrhalf.wts(:));
+
+%% polynomial interpolation
+end1 = 1;
+end2 = 1.1;
+n = 40;
+gridshifted = linspace(end1,end2,n);
+tp = [];
+tp.r = [zeros(1,n);gridshifted];
+data = ftrue(tp).'; % size = (n,1)
+fpoly = chebfun(data,'equi');
+fprime = diff(fpoly)/(end2-end1)*2;
+fprime(-1)
+
+
+%% plot
+tileplot = tiledlayout(1,3,'TileSpacing','compact');
+plotdata1 = nan(size(xx));
+plotdata1(in) = err;
+plotdata2 = nan(size(xx));
+plotdata2(in) = dudntrue;
+plotdata3 = nan(size(xx));
+plotdata3(in) = errdudn;%upcompin./dudntrue.*t.r(2,:).';
+
+
+ax1 = nexttile;
+plot(chnkrhalf); hold on; plot(chnkrhalf,'bo'); quiver(chnkrhalf,'r');
+h = pcolor(xx,yy,reshape(plotdata1,size(xx))); set(h,'EdgeColor','none');
+title('err in u'); 
+clim([-10,-1]);
+%clim([min(utrue),max(utrue)]);
+colormap(ax1,jet(100)); colorbar; axis square;
+ax2 = nexttile;
+plot(chnkrhalf); hold on; plot(chnkrhalf,'bo'); quiver(chnkrhalf,'r');
+h = pcolor(xx,yy,reshape(plotdata2,size(xx))); set(h,'EdgeColor','none');
+title('dudntrue'); clim([min(utrue),max(utrue)]); 
+clim([-18,-6]);
+colormap(ax2,jet(100)); colorbar; axis square;
+ax3 = nexttile;
+plot(chnkrhalf); hold on; plot(chnkrhalf,'bo'); quiver(chnkrhalf,'r');
+h = pcolor(xx,yy,reshape(plotdata3,size(xx))); set(h,'EdgeColor','none');
+clim([-10,-1]);
+%clim([1,3]);
+colormap(ax3,jet(100)); colorbar; axis square;
+title('upcompin');
+
+title(tileplot,'Laplace Interior Dirichlet in n-dimension');
+
diff --git a/chunkie/demo_axissymlap_analytic.m b/chunkie/demo_axissymlap_analytic.m
new file mode 100644
index 00000000..52c7de18
--- /dev/null
+++ b/chunkie/demo_axissymlap_analytic.m
@@ -0,0 +1,128 @@
+clearvars; clc;
+iftorus = 0;
+cparams = [];
+cparams.eps = 1.0e-10;
+cparams.nover = 1;
+if ~iftorus % sphere
+    cparams.ta = -pi/2;
+    cparams.tb = pi/2;
+    center = [0;2];
+    cparams.ifclosed = false;
+else % torus
+    cparams.ta = 0;
+    cparams.tb = 2*pi;
+    center = [3;0];
+    cparams.ifclosed = true;
+end
+cparams.maxchunklen = 0.5;
+radius = 1;
+fcurve = @(t) radius*[cos(t(:).'); sin(t(:).')];
+
+chnkrhalf = chunkerfunc(fcurve,cparams);
+chnkrhalf = chnkrhalf.move(-center); 
+
+figure(1);clf;plot(chnkrhalf);hold on;plot(chnkrhalf,'bo');quiver(chnkrhalf,'y');
+
+% analytic solution in n-dimension
+ndim = 5;
+ftrue = @(s) s.r(1,:).^2-(ndim-1)*s.r(2,:).^2; % analytical solution
+fdn = @(s) -2*(ndim-1)*s.r(2,:); % (2*x1,2*x2,...,-2*(ndim-1)*z)*(0,0,1) = -2*(ndim-1)*z
+fdnn = @(s) -2*(ndim-1)*ones(size(s.r(2,:)));
+
+rhs = 2*ftrue(chnkrhalf).';
+%kerns = kernel('axissymlaplace','s',ndim);
+kerndp = kernel('axissymlaplace','dp',ndim);
+kernd = kernel('axissymlaplace','d',ndim);
+D1 = 2*kernd;
+mat = chunkermat(chnkrhalf,D1);
+A = mat+eye(chnkrhalf.npt);
+sig = gmres(A, rhs, [], 1e-14, 200);
+% place targets
+x1 = linspace(center(1),center(1)+radius,200);
+x2 = linspace(center(2)-radius,center(2)+radius,200);
+[xx,yy] = meshgrid(x1,x2);
+t = [];
+t.r = [xx(:).'; yy(:).'];
+in = chunkerinterior(chnkrhalf,{x1,x2});
+t.r = t.r(:,in(:));
+ucomp = kernd.eval(chnkrhalf,t)*(sig.*chnkrhalf.wts(:));
+utrue = ftrue(t).';
+err = log10(abs((ucomp-utrue)./max(abs(utrue))));
+
+%% calculate derivative via finite difference
+eps = 1e-8;
+te = [];
+te.r = [xx(:).'; yy(:).']+[0;1]*eps;
+te.r = te.r(:,in(:));
+te.n = zeros(size(t.r))+[0;1];
+ucompe = kerndp.eval(chnkrhalf,te)*(sig.*chnkrhalf.wts(:));
+dudn_fd = (ucompe-ucomp)./eps;
+%% true derivative
+dudntrue = fdn(t).';
+errdudn_fd = log10(abs((dudntrue-dudn_fd)./dudntrue));
+
+%% calculate derivative via Dp inside domain
+t.n = zeros(size(t.r))+[0;1];
+upcompin = kerndp.eval(chnkrhalf,t)*(sig.*chnkrhalf.wts(:));
+errdudn = log10(abs((dudntrue-upcompin)./dudntrue));
+
+%% calculate derivative via Dp on boundary
+%uptrue = [2,-4]*(chnkrhalf.r(:,:).*chnkrhalf.n(:,:));
+%targpt = []; targpt.r = [0;1]; targpt.n = [0;1];
+%opts = []; opts.sing = 'hs';
+%mat = chunkermat(chnkrhalf,Dp,opts);
+%upcomp = mat*(sig.*chnkrhalf.wts(:));
+
+%% polynomial interpolation
+end1 = 1;
+end2 = 1.1;
+n = 40;
+gridshifted = linspace(end1,end2,n);
+tp = [];
+tp.r = [zeros(1,n);gridshifted];
+data = ftrue(tp).'; % size = (n,1)
+fpoly = chebfun(data,'equi');
+fprime = diff(fpoly)/(end2-end1)*2;
+fprime(-1)
+
+
+%% plot
+tileplot = tiledlayout(1,3,'TileSpacing','compact');
+plotdata1 = nan(size(xx));
+plotdata1(in) = err;
+plotdata2 = nan(size(xx));
+plotdata2(in) = dudntrue;
+plotdata3 = nan(size(xx));
+plotdata3(in) = errdudn;%upcompin./dudntrue.*t.r(2,:).';
+
+
+ax1 = nexttile;
+plot(chnkrhalf); hold on; plot(chnkrhalf,'bo'); quiver(chnkrhalf,'r');
+h = pcolor(xx,yy,reshape(plotdata1,size(xx))); set(h,'EdgeColor','none');
+title('err in u'); 
+clim([-10,-1]);
+%clim([min(utrue),max(utrue)]);
+colormap(ax1,jet(100)); colorbar; axis square;
+ax2 = nexttile;
+plot(chnkrhalf); hold on; plot(chnkrhalf,'bo'); quiver(chnkrhalf,'r');
+h = pcolor(xx,yy,reshape(plotdata2,size(xx))); set(h,'EdgeColor','none');
+title('dudntrue'); clim([min(utrue),max(utrue)]); 
+clim([-18,-6]);
+colormap(ax2,jet(100)); colorbar; axis square;
+ax3 = nexttile;
+plot(chnkrhalf); hold on; plot(chnkrhalf,'bo'); quiver(chnkrhalf,'r');
+h = pcolor(xx,yy,reshape(plotdata3,size(xx))); set(h,'EdgeColor','none');
+clim([-10,-1]);
+%clim([1,3]);
+colormap(ax3,jet(100)); colorbar; axis square;
+title('upcompin');
+
+title(tileplot,'Laplace Interior Dirichlet in n-dimension');
+
+%{
+figure(2); clf;
+zval = chnkrhalf.r(2,:).';
+plot(zval,upcomp,'DisplayName','upcomp'); hold on;
+plot(zval,uptrue,'DisplayName','uptrue'); legend;
+title('Dp(targ) for targ on boundary against z(targ)');
+%}
\ No newline at end of file