[t,u,pd,nrowpd,rpar,ipar]=djac(t,u,pd,nrowpd,rpar,ipar);

function [t,u,pd,nrowpd,rpar,ipar]=djac(t,u,pd,nrowpd,rpar,ipar);
persistent r r5 rsq u1sq u1u2 u2sq ; 

if isempty(r), r=0; end;
if isempty(r5), r5=0; end;
if isempty(rsq), rsq=0; end;
if isempty(u1sq), u1sq=0; end;
if isempty(u2sq), u2sq=0; end;
if isempty(u1u2), u1u2=0; end;
u_shape=size(u);u=reshape(u,1,[]);
pd_shape=size(pd);pd=reshape([pd(:).',zeros(1,ceil(numel(pd)./prod([nrowpd])).*prod([nrowpd])-numel(pd))],nrowpd,[]);
rpar_shape=size(rpar);rpar=reshape(rpar,1,[]);
ipar_shape=size(ipar);ipar=reshape(ipar,1,[]);
u1sq = u(1).*u(1);
u2sq = u(2).*u(2);
u1u2 = u(1).*u(2);
rsq = u1sq + u2sq;
r = sqrt(rsq);
r5 = rsq.*rsq.*r;
pd(3,1) =(3.0d0.*u1sq-rsq)./r5;
pd(4,1) = 3.0d0.*u1u2./r5;
pd(3,2) = pd(4,1);
pd(4,2) =(3.0d0.*u2sq-rsq)./r5;
pd(1,3) = 1.0d0;
pd(2,4) = 1.0d0;
u_shape=zeros(u_shape);u_shape(:)=u(1:numel(u_shape));u=u_shape;
pd_shape=zeros(pd_shape);pd_shape(:)=pd(1:numel(pd_shape));pd=pd_shape;
rpar_shape=zeros(rpar_shape);rpar_shape(:)=rpar(1:numel(rpar_shape));rpar=rpar_shape;
ipar_shape=zeros(ipar_shape);ipar_shape(:)=ipar(1:numel(ipar_shape));ipar=ipar_shape;
end %subroutine djac