function [a,ad,al,ar,all,arr]=stiff2(ni,x,p,d,Dd,c,r,k,bctype)
disp('stiff')
h=x(2:ni+1)-x(1:ni);
hcopy=h;
dcopy=d;
Ddcopy=Dd;
ccopy=c;
rcopy=r;
%
%setting constant in streamline diffusion test quantity: v+chv'
%
csd=zeros(size(c));
if p==2; hh(1:2:2*ni-1)=h; hh(2:2:2*ni)=h; end
if max(hh-d./(abs(c)+eps*ones(size(c))))>0;
csd=c./(abs(c)+eps*ones(size(c)))/2;
disp('using streamline diffusion modification');
end
%
if p==1;
%
ad(1)=d(1)/h(1)+c(1)/2+r(1)*h(1)/3;
ad(1)=ad(1)+csd(1)*(-Dd(1)+c(1)+r(1)*h(1)/2);
for i=2:ni-1;
ad(i-1)=ad(i-1)+d(i)/h(i)-c(i)/2+r(i)*h(i)/3;
ad(i-1)=ad(i-1)+csd(i)*(-Dd(i)+c(i)-r(i)*h(i)/2);
al(i-1)=-d(i)/h(i)-c(i)/2+r(i)*h(i)/6;
al(i-1)=al(i-1)+csd(i)*(Dd(i)-c(i)+r(i)*h(i)/2);
ar(i-1)=-d(i)/h(i)+c(i)/2+r(i)*h(i)/6;
ar(i-1)=ar(i-1)+csd(i)*(Dd(i)-c(i)-r(i)*h(i)/2);
ad(i)=d(i)/h(i)+c(i)/2+r(i)*h(i)/3;
ad(i)=ad(i)+csd(i)*(-Dd(i)+c(i)+r(i)*h(i)/2);
end
ad(ni-1)=ad(ni-1)+d(ni)/h(ni)-c(ni)/2+r(ni)*h(ni)/3;
ad(ni-1)=ad(ni-1)+csd(ni)*(-Dd(ni)+c(ni)-r(ni)*h(ni)/2);
%
if bctype(1)==1
ad=[d(1)/h(1)-c(1)/2+r(1)*h(1)/3+k(1)+csd(1)*(-Dd(1)+c(1)-r(1)*h(1)/2),ad];
al=[-d(1)/h(1)-c(1)/2+r(1)*h(1)/6+csd(1)*(Dd(1)-c(1)+r(1)*h(1)/2),al];
ar=[-d(1)/h(1)+c(1)/2+r(1)*h(1)/6+csd(1)*(Dd(1)-c(1)-r(1)*h(1)/2),ar];
end
if bctype(2)==1
ad=[ad,d(ni)/h(ni)+c(ni)/2+r(ni)*h(ni)/3+k(2)+csd(ni)*(-Dd(ni)+c(ni)+r(ni)*h(ni)/2)];
al=[al,-d(ni)/h(ni)-c(ni)/2+r(ni)*h(ni)/6+csd(ni)*(Dd(ni)-c(ni)+r(ni)*h(ni)/2)];
end
%
elseif p==2;
%
gauss1=1/2-1/sqrt(3)/2;
gauss2=1/2+1/sqrt(3)/2;
v1=4*gauss1*(1-gauss1);
v2=4*gauss2*(1-gauss2);
dv1=4-8*gauss1;
dv2=4-8*gauss2;
wl1=(1-gauss2)*(1-2*gauss2);
wl2=(1-gauss1)*(1-2*gauss1);
dwl1=-4*gauss2+3;
dwl2=-4*gauss1+3;
wr1=(1-gauss1)*(1-2*gauss1);
wr2=(1-gauss2)*(1-2*gauss2);
dwr1=4*gauss1-3;
dwr2=4*gauss2-3;
%
%first interval
%
h=hcopy(1);
d1=dcopy(1);
d2=dcopy(2);
Dd1=Ddcopy(1);
Dd2=Ddcopy(2);
c1=ccopy(1);
c2=ccopy(2);
r1=rcopy(1);
r2=rcopy(2);
csd1=csd(1);
csd2=csd(2);
ad(1)=      (d1*dv1 *dv1 /h+c1*dv1 *v1 +r1*v1 *v1 *h)/2;
ad(1)=ad(1)+(d2*dv2 *dv2 /h+c2*dv2 *v2 +r2*v2 *v2 *h)/2;
ad(1)=ad(1)+csd1*(d1*8/h-Dd1*dv1+c1*dv1+r1*v1*h)*dv1/2;
ad(1)=ad(1)+csd2*(d2*8/h-Dd2*dv2+c2*dv2+r2*v2*h)*dv2/2;
ad(2)=      (d1*dwl1*dwl1/h+c1*dwl1*wl1+r1*wl1*wl1*h)/2;
ad(2)=ad(2)+(d2*dwl2*dwl2/h+c2*dwl2*wl2+r2*wl2*wl2*h)/2;
ad(2)=ad(2)+csd1*(-d1*4/h-Dd1*dwl1+c1*dwl1+r1*wl1*h)*dwl1/2;
ad(2)=ad(2)+csd2*(-d2*4/h-Dd2*dwl2+c2*dwl2+r2*wl2*h)*dwl2/2;
al(1)=      (d1*dv1 *dwl1/h+c1*dv1 *wl1+r1*v1 *wl1*h)/2;
al(1)=al(1)+(d2*dv2 *dwl2/h+c2*dv2 *wl2+r2*v2 *wl2*h)/2;
al(1)=al(1)+csd1*(d1*8/h-Dd1*dv1+c1*dv1+r1*v1*h)*dwl1/2;
al(1)=al(1)+csd2*(d2*8/h-Dd2*dv2+c2*dv2+r2*v2*h)*dwl2/2;
ar(1)=      (d1*dwl1*dv1 /h+c1*dwl1*v1 +r1*wl1*v1 *h)/2;
ar(1)=ar(1)+(d2*dwl2*dv2 /h+c2*dwl2*v2 +r2*wl2*v2 *h)/2;
ar(1)=ar(1)+csd1*(-d1*4/h-Dd1*dwl1+c1*dwl1+r1*wl1*h)*dv1/2;
ar(1)=ar(1)+csd2*(-d2*4/h-Dd2*dwl2+c2*dwl2+r2*wl2*h)*dv2/2;
%
%interior elements
%
for i=2:ni-1;
h=hcopy(i);
d1=dcopy(2*i-1);
d2=dcopy(2*i);
Dd1=Ddcopy(2*i-1);
Dd2=Ddcopy(2*i);
c1=ccopy(2*i-1);
c2=ccopy(2*i);
r1=rcopy(2*i-1);
r2=rcopy(2*i);
csd1=csd(2*i-1);
csd2=csd(2*i);
ad(2*i-2)=ad(2*i-2)  +(d1*dwr1*dwr1/h+c1*dwr1*wr1+r1*wr1*wr1*h)/2;
ad(2*i-2)=ad(2*i-2)  +(d2*dwr2*dwr2/h+c2*dwr2*wr2+r2*wr2*wr2*h)/2;
ad(2*i-2)=ad(2*i-2)+csd1*(-d1*4/h-Dd1*dwr1+c1*dwr1+r1*wr1*h)*dwr1/2;
ad(2*i-2)=ad(2*i-2)+csd2*(-d2*4/h-Dd2*dwr2+c2*dwr2+r2*wr2*h)*dwr2/2;
ad(2*i-1)=            (d1*dv1 *dv1 /h+c1*dv1 *v1 +r1*v1 *v1 *h)/2;
ad(2*i-1)=ad(2*i-1)  +(d2*dv2 *dv2 /h+c2*dv2 *v2 +r2*v2 *v2 *h)/2;
ad(2*i-1)=ad(2*i-1)+csd1*(d1*8/h-Dd1*dv1+c1*dv1+r1*v1*h)*dv1/2;
ad(2*i-1)=ad(2*i-1)+csd2*(d2*8/h-Dd2*dv2+c2*dv2+r2*v2*h)*dv2/2;
ad(2*i)=              (d1*dwl1*dwl1/h+c1*dwl1*wl1+r1*wl1*wl1*h)/2;
ad(2*i)=ad(2*i)      +(d2*dwl2*dwl2/h+c2*dwl2*wl2+r2*wl2*wl2*h)/2;
ad(2*i)=ad(2*i)+csd1*(-d1*4/h-Dd1*dwl1+c1*dwl1+r1*wl1*h)*dwl1/2;
ad(2*i)=ad(2*i)+csd2*(-d2*4/h-Dd2*dwl2+c2*dwl2+r2*wl2*h)*dwl2/2;
al(2*i-2)=            (d1*dwr1*dv1 /h+c1*dwr1*v1 +r1*wr1*v1 *h)/2;
al(2*i-2)=al(2*i-2)  +(d2*dwr2*dv2 /h+c2*dwr2*v2 +r2*wr2*v2 *h)/2;
al(2*i-2)=al(2*i-2)+csd1*(-d1*4/h-Dd1*dwr1+c1*dwr1+r1*wr1*h)*dv1/2;
al(2*i-2)=al(2*i-2)+csd2*(-d2*4/h-Dd2*dwr2+c2*dwr2+r2*wr2*h)*dv2/2;
ar(2*i-2)=            (d1*dv1 *dwr1/h+c1*dv1*wr1 +r1*v1 *wr1*h)/2;
ar(2*i-2)=ar(2*i-2)  +(d2*dv2 *dwr2/h+c2*dv2*wr2 +r2*v2 *wr2*h)/2;
ar(2*i-2)=ar(2*i-2)+csd1*(d1*8/h-Dd1*dv1+c1*dv1+r1*v1*h)*dwr1/2;
ar(2*i-2)=ar(2*i-2)+csd2*(d2*8/h-Dd2*dv2+c2*dv2+r2*v2*h)*dwr2/2;
al(2*i-1)=            (d1*dv1*dwl1 /h+c1*dv1*wl1 +r1*v1 *wl1*h)/2;
al(2*i-1)=al(2*i-1)  +(d2*dv2*dwl2 /h+c2*dv2 *wl2+r2*v2 *wl2*h)/2;
al(2*i-1)=al(2*i-1)+csd1*(d1*8/h-Dd1*dv1+c1*dv1+r1*v1*h)*dwl1/2;
al(2*i-1)=al(2*i-1)+csd2*(d2*8/h-Dd2*dv2+c2*dv2+r2*v2*h)*dwl2/2;
ar(2*i-1)=            (d1*dwl1*dv1 /h+c1*dwl1*v1 +r1*wl1*v1 *h)/2;
ar(2*i-1)=ar(2*i-1)  +(d2*dwl2*dv2 /h+c2*dwl2*v2 +r2*wl2*v2 *h)/2;
ar(2*i-1)=ar(2*i-1)+csd1*(-d1*4/h-Dd1*dwl1+c1*dwl1+r1*wl1*h)*dv1/2;
ar(2*i-1)=ar(2*i-1)+csd2*(-d2*4/h-Dd2*dwl2+c2*dwl2+r2*wl2*h)*dv2/2;
all(2*i-2)=           (d1*dwr1*dwl1/h+c1*dwr1*wl1+r1*wr1*wl1*h)/2;
all(2*i-2)=all(2*i-2)+(d2*dwr2*dwl2/h+c2*dwr2*wl2+r2*wr2*wl2*h)/2;
all(2*i-2)=all(2*i-2)+csd1*(-d1*4/h-Dd1*dwr1+c1*dwr1+r1*wr1*h)*dwl1/2;
all(2*i-2)=all(2*i-2)+csd2*(-d2*4/h-Dd2*dwr2+c2*dwr2+r2*wr2*h)*dwl2/2;
arr(2*i-2)=           (d1*dwl1*dwr1/h+c1*dwl1*wr1+r1*wl1*wr1*h)/2;
arr(2*i-2)=arr(2*i-2)+(d2*dwl2*dwr2/h+c2*dwl2*wr2+r2*wl2*wr2*h)/2;
arr(2*i-2)=arr(2*i-2)+csd1*(-d1*4/h-Dd1*dwl1+c1*dwl1+r1*wl1*h)*dwr1/2;
arr(2*i-2)=arr(2*i-2)+csd2*(-d2*4/h-Dd2*dwl2+c2*dwl2+r2*wl2*h)*dwr2/2;
end
%
%last interval
%
h=hcopy(ni);
d1=dcopy(2*ni-1);
d2=dcopy(2*ni);
Dd1=Ddcopy(2*ni-1);
Dd2=Ddcopy(2*ni);
c1=ccopy(2*ni-1);
c2=ccopy(2*ni);
r1=rcopy(2*ni-1);
r2=rcopy(2*ni);
csd1=csd(2*ni-1);
csd2=csd(2*ni);
ad(2*ni-2)=ad(2*ni-2)+(d1*dwr1*dwr1/h+c1*dwr1*wr1+r1*wr1*wr1*h)/2;
ad(2*ni-2)=ad(2*ni-2)+(d2*dwr2*dwr2/h+c2*dwr2*wr2+r2*wr2*wr2*h)/2;
ad(2*ni-2)=ad(2*ni-2)+csd1*(-d1*4/h-Dd1*dwr1+c1*dwr1+r1*wr1*h)*dwr1/2;
ad(2*ni-2)=ad(2*ni-2)+csd2*(-d2*4/h-Dd2*dwr2+c2*dwr2+r2*wr2*h)*dwr2/2;
ad(2*ni-1)=           (d1*dv1 *dv1 /h+c1*dv1 *v1 +r1*v1 *v1 *h)/2;
ad(2*ni-1)=ad(2*ni-1)+(d2*dv2 *dv2 /h+c2*dv2 *v2 +r2*v2 *v2 *h)/2;
ad(2*ni-1)=ad(2*ni-1)+csd1*(d1*8/h-Dd1*dv1+c1*dv1+r1*v1*h)*dv1/2;
ad(2*ni-1)=ad(2*ni-1)+csd2*(d2*8/h-Dd2*dv2+c2*dv2+r2*v2*h)*dv2/2;
al(2*ni-2)=           (d1*dwr1*dv1 /h+c1*dwr1*v1 +r1*wr1*v1 *h)/2;
al(2*ni-2)=al(2*ni-2)+(d2*dwr2*dv2 /h+c2*dwr2*v2 +r2*wr2*v2 *h)/2;
al(2*ni-2)=al(2*ni-2)+csd1*(-d1*4/h-Dd1*dwr1+c1*dwr1+r1*wr1*h)*dv1/2;
al(2*ni-2)=al(2*ni-2)+csd2*(-d2*4/h-Dd2*dwr2+c2*dwr2+r2*wr2*h)*dv2/2;
ar(2*ni-2)=           (d1*dv1*dwr1 /h+c1*dv1 *wr1+r1*v1 *wr1*h)/2;
ar(2*ni-2)=ar(2*ni-2)+(d2*dv2*dwr2 /h+c2*dv2 *wr2+r2*v2 *wr2*h)/2;
ar(2*ni-2)=ar(2*ni-2)+csd1*(d1*8/h-Dd1*dv1+c1*dv1+r1*v1*h)*dwr1/2;
ar(2*ni-2)=ar(2*ni-2)+csd2*(d2*8/h-Dd2*dv2+c2*dv2+r2*v2*h)*dwr2/2;
all=[all,0];
arr=[arr,0];
%
if bctype(1)==1
h=hcopy(1);
d1=d(1);
d2=d(2);
c1=c(1);
c2=c(2);
r1=r(1);
r2=r(2);
Dd1=Dd(1);
Dd2=Dd(2);
csd1=csd(1);
csd2=csd(2);
ad=[0,ad];
al=[0,al];
ar=[0,ar];
all=[0,all];
arr=[0,arr];
ad(1)=(d1*dwr1*dwr1/h+c1*dwr1*wr1+r1*wr1*wr1*h)/2;
ad(1)=ad(1)+(d2*dwr2*dwr2/h+c2*dwr2*wr2+r2*wr2*wr2*h)/2;
ad(1)=ad(1)+csd1*(-d1*4/h-Dd1*dwr1+c1*dwr1+r1*wr1*h)*dwr1/2;
ad(1)=ad(1)+csd2*(-d2*4/h-Dd2*dwr2+c2*dwr2+r2*wr2*h)*dwr2/2;
al(1)=(d1*dwr1*dv1/h+c1*dwr1*v1+r1*wr1*v1*h)/2;
al(1)=al(1)+(d2*dwr2*dv2/h+c2*dwr2*v2+r2*wr2*v2*h)/2;
al(1)=al(1)+csd1*(-d1*4/h-Dd1*dwr1+c1*dwr1+r1*wr1*h)*dv1/2;
al(1)=al(1)+csd2*(-d2*4/h-Dd2*dwr2+c2*dwr2+r2*wr2*h)*dv2/2;
ar(1)=(d1*dv1*dwr1/h+c1*dv1*wr1+r1*v1*wr1*h)/2;
ar(1)=ar(1)+(d2*dv2*dwr2/h+c2*dv2*wr2+r2*v2*wr2*h)/2;
ar(1)=ar(1)+csd1*(d1*8/h-Dd1*dv1+c1*dv1+r1*v1*h)*dwr1/2;
ar(1)=ar(1)+csd2*(d2*8/h-Dd2*dv2+c2*dv2+r2*v2*h)*dwr2/2;
arr(1)=(d1*dwl1*dwr1/h+c1*dwl1*wr1+r1*wl1*wr1*h)/2;
arr(1)=arr(1)+(d2*dwl2*dwr2/h+c2*dwl2*wr2+r2*wl2*wr2*h)/2;
arr(1)=arr(1)+csd1*(-d1*4/h-Dd1*dwl1+c1*dwl1+r1*wl1*h)*dwr1/2;
arr(1)=arr(1)+csd2*(-d2*4/h-Dd2*dwl2+c2*dwl2+r2*wl2*h)*dwr2/2;
all(1)=(d1*dwr1*dwl1/h+c1*dwr1*wl1+r1*wr1*wl1*h)/2;
all(1)=all(1)+(d2*dwr2*dwl2/h+c2*dwr2*wl2+r2*wr2*wl2*h)/2;
all(1)=all(1)+csd1*(-d1*4/h-Dd1*dwr1+c1*dwr1+r1*wr1*h)*dwl1/2;
all(1)=all(1)+csd2*(-d2*4/h-Dd2*dwr2+c2*dwr2+r2*wr2*h)*dwl2/2;
end
if bctype(2)==1
h=hcopy(ni);
d1=d(2*ni-1);
d2=d(2*ni);
Dd1=Dd(2*ni-1);
Dd2=Dd(2*ni);
c1=c(2*ni-1);
c2=c(2*ni);
r1=r(2*ni-1);
r2=r(2*ni);
csd1=csd(2*ni-1);
csd2=csd(2*ni);
ad=[ad,0];
al=[al,0];
ar=[ar,0];
l=length(ad);
ad(l)=(d1*dwl1*dwl1/h+c1*dwl1*wl1+r1*wl1*wl1*h)/2;
ad(l)=ad(l)+(d2*dwl2*dwl2/h+c2*dwl2*wl2+r2*wl2*wl2*h)/2;
ad(l)=ad(l)+csd1*(-d1*4/h-Dd1*dwl1+c1*dwl1+r1*wl1*h)*dwl1/2;
ad(l)=ad(l)+csd2*(-d2*4/h-Dd2*dwl2+c2*dwl2+r2*wl2*h)*dwl2/2;
al(l-1)=(d1*dwr1*dv1/h+c1*dwr1*v1+r1*wr1*v1*h)/2;
al(l-1)=al(l-1)+(d2*dwr2*dv2/h+c2*dwr2*v2+r2*wr2*v2*h)/2;
al(l-1)=al(l-1)+csd1*(d1*8/h-Dd1*dv1+c1*dv1+r1*v1*h)*dwl1/2;
al(l-1)=al(l-1)+csd2*(d2*8/h-Dd2*dv2+c2*dv2+r2*v2*h)*dwl2/2;
ar(l-1)=(d1*dwl1*dv1/h+c1*dwl1*v1+r1*wl1*v1*h)/2;
ar(l-1)=ar(l-1)+(d2*dwl2*dv2/h+c2*dwl2*v2+r2*wl2*v2*h)/2;
ar(l-1)=ar(l-1)+csd1*(-d1*4/h-Dd1*dwl1+c1*dwl1+r1*wl1*h)*dv1/2;
ar(l-1)=ar(l-1)+csd2*(-d2*4/h-Dd2*dwl2+c2*dwl2+r2*wl2*h)*dv2/2;
all(l-2)=(d1*dwr1*dwl1/h+c1*dwr1*wl1+r1*wr1*wl1*h)/2;
all(l-2)=all(l-2)+(d2*dwr2*dwl2/h+c2*dwr2*wl2+r2*wr2*wl2*h)/2;
all(l-2)=all(l-2)+csd1*(-d1*4/h-Dd1*dwr1+c1*dwr1+r1*wr1*h)*dwl1/2;
all(l-2)=all(l-2)+csd2*(-d2*4/h-Dd2*dwr2+c2*dwr2+r2*wr2*h)*dwl2/2;
arr(l-2)=(d1*dwl1*dwr1/h+c1*dwl1*wr1+r1*wl1*wr1*h)/2;
arr(l-2)=arr(l-2)+(d2*dwl2*dwl2/h+c2*dwl2*wr2+r2*wl2*wr2*h)/2;
arr(l-2)=arr(l-2)+csd1*(-d1*4/h-Dd1*dwl1+c1*dwl1+r1*wl1*h)*dwr1/2;
arr(l-2)=arr(l-2)+csd2*(-d2*4/h-Dd2*dwl2+c2*dwl2+r2*wl2*h)*dwr2/2;
end
ar=[0 ar];
arr=[0 0 arr];
al=[al 0];
all=[all 0 0];
l=length(ad);
a=full(spdiags([all' al' ad' ar' arr'],[-2 -1 0 1 2],l,l))
end %of if p branch
