function [z,dz,d2z,d3z,dpz,dxpz,d2pz]=dmult(v,dv,d2v,d3v,dpv,dxpv,d2pv, ...
                                            w,dw,d2w,d3w,dpw,dxpw,d2pw);
%DMULT   Calculates 3rd-order derivatives for v*w when v,w's are known.
%        [z,dz,d2z,d3z,dpz,dxpz,d2pz]=dmult(v,dv,d2v,d3v,dpv,dxpv,d2pv, ...
%        w,dw,d2w,d3w,dpw,dxpw,d2pw) calculates derivatives dz/dx, 
%        d^2z/dx^2,d^3z/dx^3,dz/dp,d^2z/dxdp,d^3z/dx^2dp for z=v(x,p)*w(x,p).

%        Ellen McGrattan, 5-20-91
%        Revised, 5-20-91, ERM

e   = eye(size(d2v));
z   = v*w;
dz  = v*dw+w*dv;
d2z = v*d2w+dv*dw'+dw*dv'+w*d2v;
d3z = v*d3w+vec(d2w)*dv'+(kron(e,dv)+kron(dv,e))*d2w+ ...
      (kron(e,dw)+kron(dw,e))*d2v+w*d3v+vec(d2v)*dw';
dpz = v*dpw+w*dpv;
dxpz= v*dxpw+dv*dpw'+w*dxpv+dw*dpv';
d2pz= v*d2pw+vec(d2w)*dpv'+(kron(e,dv)+kron(dv,e))*dxpw+ ...
      (kron(e,dw)+kron(dw,e))*dxpv+w*d2pv+vec(d2v)*dpw';
