[t,y,yprime,pd,cj,rpar,ipar]=ddjac2(t,y,yprime,pd,cj,rpar,ipar);

function [t,y,yprime,pd,cj,rpar,ipar]=ddjac2(t,y,yprime,pd,cj,rpar,ipar);
persistent alph1 alph2 firstCall j mband ml mu neq ng ; if isempty(firstCall),firstCall=1;end; 

ipar_shape=size(ipar);ipar=reshape(ipar,1,[]);
y_shape=size(y);y=reshape(y,1,[]);
yprime_shape=size(yprime);yprime=reshape(yprime,1,[]);
pd_orig=pd;pd_shape=[11,25];pd=reshape([pd_orig(1:min(prod(pd_shape),numel(pd_orig))),zeros(1,max(0,prod(pd_shape)-numel(pd_orig)))],pd_shape);
rpar_shape=size(rpar);rpar=reshape(rpar,1,[]);
if isempty(j), j=0; end;
if isempty(mband), mband=0; end;
if isempty(ml), ml=0; end;
if isempty(mu), mu=0; end;
if isempty(neq), neq=0; end;
if isempty(ng), ng=0; end;
if isempty(alph1), alph1=0; end;
if isempty(alph2), alph2=0; end;
if firstCall,   alph1=[1.0d0];  end;
if firstCall,  alph2=[1.0d0];  end;
if firstCall,  ng=[5];  end;
if firstCall,   ml=[5];  end;
if firstCall,  mu=[0];  end;
if firstCall,  neq=[25];  end;
firstCall=0;
mband = fix(ml + mu + 1);
for j = 1 : neq;
pd(mband,j) = -2.0d0 - cj;
pd(mband+1,j) = alph1;
pd(mband+2,j) = 0.0d0;
pd(mband+3,j) = 0.0d0;
pd(mband+4,j) = 0.0d0;
pd(mband+5,j) = alph2;
end; j = fix(neq+1);
for j = 1 : ng: neq ;
pd(mband+1,j) = 0.0d0;
end; j = fix(neq +1);
ipar_shape=zeros(ipar_shape);ipar_shape(:)=ipar(1:numel(ipar_shape));ipar=ipar_shape;
y_shape=zeros(y_shape);y_shape(:)=y(1:numel(y_shape));y=y_shape;
yprime_shape=zeros(yprime_shape);yprime_shape(:)=yprime(1:numel(yprime_shape));yprime=yprime_shape;
pd_orig(1:prod(pd_shape))=pd;pd=pd_orig;
rpar_shape=zeros(rpar_shape);rpar_shape(:)=rpar(1:numel(rpar_shape));rpar=rpar_shape;
end %subroutine ddjac2