function [a,b]=stiffload(ni,x,p,d,Dd,c,r,f,k,g,bctype)
disp('stiffload')
h=x(2:ni+1)-x(1:ni);
%
if p==1;
%
csd=zeros(h);
if max(h-d./(abs(c)+eps*ones(c)))>0;
csd=c./(abs(c)+eps*ones(c))/2;
disp('using streamline diffusion');
end
%
wr=0.5*ones(h);
wl=0.5*ones(h);
dwr=-ones(h)./h;
dwl=ones(h)./h;
%
wrwl=(d.*dwr.*dwl+c.*dwr.*wl+r.*wr.*wl+(-Dd.*dwr+c.*dwr+r.*wr).*csd.*h.*dwl).*h;
wlwr=(d.*dwl.*dwr+c.*dwl.*wr+r.*wl.*wr+(-Dd.*dwl+c.*dwl+r.*wl).*csd.*h.*dwr).*h;
wrwr=(d.*dwr.*dwr+c.*dwr.*wr+r.*wr.*wr+(-Dd.*dwr+c.*dwr+r.*wr).*csd.*h.*dwr).*h;
wlwl=(d.*dwl.*dwl+c.*dwl.*wl+r.*wl.*wl+(-Dd.*dwl+c.*dwl+r.*wl).*csd.*h.*dwl).*h;
fwr=f.*(wr+csd.*h.*dwr).*h;
fwl=f.*(wl+csd.*h.*dwl).*h;
%
ad(2:ni+1)=wlwl;
ad(1:ni)=ad(1:ni)+wrwr;
al(1:ni)=wrwl;
ar(1:ni)=wlwr;
b(2:ni+1)=fwl;
b(1:ni)=b(1:ni)+fwr;
%
ar=[ar 0 ];
al=[al 0];
for i=1:ni+1;
a(i,i)=ad(i);
a(i+1,i)=al(i);
a(i,i+1)=ar(i);
end
a=a(1:ni+1,1:ni+1);
%
elseif p==2;
%
hh(2:2:2*ni)=h; hh(1:2:2*ni-1)=h;
%
csd=zeros(hh);
if max(hh-d./(abs(c)+eps*ones(c)))>0;
csd=c./(abs(c)+eps*ones(c))/2;
disp('using streamline diffusion');
end
%
gauss=1/sqrt(3);
gauss1=(1-gauss)/2;
gauss2=(1+gauss)/2;
%
v=4*gauss1*(1-gauss1)*ones(d);
dv(1:2:2*ni-1)=(4-8*gauss1)*ones(h)./h;
dv(2:2:2*ni)=-dv(1:2:2*ni-1);
ddv=-8*ones(d);
wr(1:2:2*ni-1)=(1-gauss1)*(1-2*gauss1)*ones(h);
wr(2:2:2*ni)=(1-gauss2)*(1-2*gauss2)*ones(h);
dwr(1:2:2*ni-1)=(4*gauss1-3)*ones(h)./h;
dwr(2:2:2*ni)=(4*gauss2-3)*ones(h)./h;
ddwr=4*ones(d);
wl(1:2:2*ni-1)=wr(2:2:2*ni);
wl(2:2:2*ni)=wr(1:2:2*ni-1);
dwl(1:2:2*ni-1)=-dwr(2:2:2*ni);
dwl(2:2:2*ni)=-dwr(1:2:2*ni-1);
ddwl=4*ones(d);
%
vv=(d.*dv.*dv+c.*dv.*v+r.*v.*v+(-d.*ddv-Dd.*dv+c.*dv+r.*v).*csd.*hh.*dv).*(hh/2);
wrv=(d.*dwr.*dv+c.*dwr.*v+r.*wr.*v+(-d.*ddwr-Dd.*dwr+c.*dwr+r.*wr).*csd.*hh.*dv).*(hh/2);
wlv=(d.*dwl.*dv+c.*dwl.*v+r.*wl.*v+(-d.*ddwl-Dd.*dwl+c.*dwl+r.*wl).*csd.*hh.*dv).*(hh/2);
vwl=(d.*dv.*dwl+c.*dv.*wl+r.*v.*wl+(-d.*ddv-Dd.*dv+c.*dv+r.*v).*csd.*hh.*dwl).*(hh/2);
vwr=(d.*dv.*dwr+c.*dv.*wr+r.*v.*wr+(-d.*ddv-Dd.*dv+c.*dv+r.*v).*csd.*hh.*dwr).*(hh/2);
wrwl=(d.*dwr.*dwl+c.*dwr.*wl+r.*wr.*wl+(-d.*ddwr-Dd.*dwr+c.*dwr+r.*wr).*csd.*hh.*dwl).*(hh/2);
wlwr=(d.*dwl.*dwr+c.*dwl.*wr+r.*wl.*wr+(-d.*ddwl-Dd.*dwl+c.*dwl+r.*wl).*csd.*hh.*dwr).*(hh/2);
wrwr=(d.*dwr.*dwr+c.*dwr.*wr+r.*wr.*wr+(-d.*ddwr-Dd.*dwr+c.*dwr+r.*wr).*csd.*hh.*dwr).*(hh/2);
wlwl=(d.*dwl.*dwl+c.*dwl.*wl+r.*wl.*wl+(-d.*ddwl-Dd.*dwl+c.*dwl+r.*wl).*csd.*hh.*dwl).*(hh/2);
fv=f.*(v+csd.*hh.*dv).*(hh/2);
fwr=f.*(wr+csd.*hh.*dwr).*(hh/2);
fwl=f.*(wl+csd.*hh.*dwl).*(hh/2);
%
vv=vv(1:2:2*ni-1)+vv(2:2:2*ni);
wrv=wrv(1:2:2*ni-1)+wrv(2:2:2*ni);
wlv=wlv(1:2:2*ni-1)+wlv(2:2:2*ni);
vwl=vwl(1:2:2*ni-1)+vwl(2:2:2*ni);
vwr=vwr(1:2:2*ni-1)+vwr(2:2:2*ni);
wrwl=wrwl(1:2:2*ni-1)+wrwl(2:2:2*ni);
wlwr=wlwr(1:2:2*ni-1)+wlwr(2:2:2*ni);
wrwr=wrwr(1:2:2*ni-1)+wrwr(2:2:2*ni);
wlwl=wlwl(1:2:2*ni-1)+wlwl(2:2:2*ni);
fv=fv(1:2:2*ni-1)+fv(2:2:2*ni);
fwr=fwr(1:2:2*ni-1)+fwr(2:2:2*ni);
fwl=fwl(1:2:2*ni-1)+fwl(2:2:2*ni);
%
ad(2:2:2*ni)=vv;
ad(3:2:2*ni+1)=wlwl;
ad(1:2:2*ni-1)=ad(1:2:2*ni-1)+wrwr;
al(1:2:2*ni-1)=wrv;
al(2:2:2*ni)=vwl;
ar(1:2:2*ni-1)=vwr;
ar(2:2:2*ni)=wlv;
all(1:2:2*ni-1)=wrwl;
arr(1:2:2*ni-1)=wlwr;
b(2:2:2*ni)=fv;
b(3:2:2*ni+1)=fwl;
b(1:2:2*ni-1)=b(1:2:2*ni-1)+fwr;
%
ar=[ar 0 ];
arr=[arr 0 0];
al=[al 0];
all=[all 0 0];
for i=1:2*ni+1;
a(i,i)=ad(i);
a(i+1,i)=al(i);
a(i+2,i)=all(i);
a(i,i+1)=ar(i);
a(i,i+2)=arr(i);
end
a=a(1:2*ni+1,1:2*ni+1);
end %if p
%
ndf=p*ni-1+sum(bctype);
if bctype(1)==0;
a=a(2:p*ni+1,2:p*ni+1);
b=b(2:p*ni+1);
b=b-g(1)*a(:,1)';
elseif bctype(1)==1;
a(1,1)=a(1,1)+k(1);
b(1)=b(1)+g(1);
end
if bctype(2)==0;
a=a(1:ndf,1:ndf);
b=b(1:ndf);
b=b-g(2)*a(:,ndf)';
elseif bctype(2)==1;
a(ndf,ndf)=a(ndf,ndf)+k(2);
b(ndf)=b(ndf)+g(2);
end
