ADFEM manual:

all the files necessary are m-files and are found
further down in this ascii-file. You have to split the file into
the m-files 
adfem.m
solve.m
corr.m
corrnew.m
loads.m
stiff.m
pred.m
test.m
parameters.m
exacterr.m
the files are separated by a line like
---------new file------
and the end of file line
is
--------------end of file-----------------

---------new file------

   ADFEM is an adaptive finite element matlab code for approximate solution 
of two-point boundary value problems of the form: Find  u = u(x)  such that, 
with  D = d/dx,

	    - D(d(x)Du) + c(x)Du + a(x)u = f(x) for xmin < x < xmax,

and such that

		     dDu + ku = g    or    u = g

at x=xmin, and

     - dDu + ku = g    or    u = g

at x=xmax. 

   Here d(x) is the diffusion coefficient (positive), c(x) represents 
convection, a(x) may be interpreted as an absorption intensity, and 
f(x) is a force or load acting on the system.

   The boundary conditions may be of Dirichlet type, with the value of  u  
given, or of Neumann/Robin type, with a combination of the normal flux  ądDu
and  u  given. In a heat conduction problem, where  u  is the temperature 
and  -dDu  is the heat flux, this corresponds to prescribing the temperature 
at the boundary, prescribing the heat flux through the boundary (with k=0), 
or relating the heat flux through the boundary to the temperature jump 
through the boundary in which case  g/k  represents the exterior temperature. 
Of course, the values of  k  and  g  need not be the same at  x=a  and  x=b.
   
   The discretization method used is the most standard finite element method 
using piecewise linear trial and test functions (with a stabilizing "stream-
line diffusion modification" in case of convection dominated problems). 

   The error may be controlled in the energy, L2 or maximum norm on a user 
specified tolerance level  tol.

   The structure of the adaptive algorithm used is as follows. Once all data 
of the problem have been given, the program predicts (in a subroutine "pred") 
a suitable initial (usually uniform) computational mesh, and computes the 
corresponding finite element approximation to the solution. This solution is 
then tested (in a subroutine "test") using a so called a posteriori estimate 
for the error (in the given norm), where the residual of the approximate solu-
tion plays an important role. If the current approximate solution is found 
to have an error which is OK, it is presented as "tested and accepted" and 
the algorithm has come to a successful stop. Otherwise, a new and improved 
(corrected) computational mesh is constructed (in a subroutine "corr") using 
the information obtained in the previous test, and the above solution and 
test procedure is repeated. If we are lucky the new current solution is now 
OK. Otherwise, the adaptive procedure just described is repeated until a 
successful stop is reached. The current (local) meshsize, solution and resi-
dual are plotted throughout the adaptive process. In the error estimate used 
to control the error a certain stability factor enters which is also dis-
played. The size of this constant in some sense reflects the "difficulty" of 
the given problem. 

   To run the program, first give the command  matlab  to start the matlab 
program. Then write  adfem  at the matlab >> prompt. Input data and await 
the desired solution!

   In case you want to re-run adfem, with essentially the same data, you can 
take a short-cut by giving the new data directly at the matlab >> prompt and
then give the command  solve.  For example, if you want to re-run with a 
finer
error tolerance you first write  tol=0.001  and then  solve,  and adfem will 
resolve the given problem with the new tolerance  tol. Other data that you 
may want to reset in this way are
xmin (=a), xmax (=b), diffusion, convection, absorption, force, 
errnorm, exactflag, uexact, 
bctype, k and g.
Here the last three are 2-vectors (1x2 matrices in matlab) defining the 
boundary condition type (=0 for Dirichlet and =1 for Neumann/Robin) and the 
corresponding k and g values at x=a and x=b. xmin and xmax defines the 
computational domain. The rest of the data are text strings and has to be 
gives within single quotes like, for example, force='1/x'. The error norm 
can be chosen by errnorm='en' (for energy), ='l2' (for L2) or ='ma' (for 
the maximum norm). Finally there is
an option where you are given the possibility to check the efficiency and 
reliability of the error control by providing the exact solution uexact. 
The program will then compute the true error (exacterror) for comparison 
with the estimated one (error). exactflag should be set to 'y' (for yes) 
in this case. 

Good luck !
-----------new file---------
clear
disp('This program solves the diffusion-convection-absorbtion equation')
disp('')
disp('    -D(dDu) + cDu + au = f      (D=d/dx)')
disp('')
disp('in the domain  [xmin,xmax],')
disp('with boundary conditions')
disp('')
disp('     dDu + ku = g   or   u = g')
disp('')
disp('at  x = xmin  and  x = xmax,  respectively,')
disp('')
disp('and with energy, L2 or maximum norm error at most tol.')
%
%data input
%
xmin=input('xmin = ');
xmax=input('xmax = ');
%
diffusion=input('d(x) = ','s');
convection=input('c(x) = ','s');
absorbtion=input('a(x) = ','s');
force=input('f(x) = ','s');
%
k=[0 0];
input('b.c. type at  x = xmin,  n (Neuman/Robin) or d (Dirichlet) ? ','s');
if ans=='d'; bctype(1)=0; else; bctype(1)=1; end
if bctype(1)==1
k(1)=input('k = ');
end
g(1)=input('g = ');
input('b.c. type at  x = xmax,  n or d ? ','s');
if ans=='d'; bctype(2)=0; else; bctype(2)=1; end
if bctype(2)==1
k(2)=input('k = ');
end
g(2)=input('g = ');
%
errnorm=input('error norm ?  en (energy), l2  or  ma (maximum) ? ','s');
if convection=='0';
%disp('hope absorbtion term  a  is not (too) negative')
else
if errnorm=='en';
disp('problem non symmetric - energy norm does not exist')
disp('choose another error norm')
errnorm=input('l2  or  ma ? ','s');
end
end
tol=input('tol = ');
%
disp('is the exact solution known?  y  or  n ? ')
exactflag=input(' ','s');
if exactflag=='y';
disp('give  u(x) ')
uexact=input('u(x) = ','s');
end
%
run=0;
solve

--------------new file ------------
clear d Dd c r f xold x xcopy ni ndf h hh ad al ar adcopy b bcopy
clear uold u res globalerr ok ans i uex Duex exacterror
if run==0
close
figure('Position',[500 40 500 850])
else
clg
end
parameters
%
[ni,x]=pred(xmin,xmax,force,errnorm,tol);
%

if exactflag=='y';
xcopy=x;
%for i=1:201;
%x=xmin+(i-1)*(xmax-xmin)/200;
%uexplot(i)=eval(uexact);
%end
xplot=[xmin:(xmax-xmin)/200:xmax];
x=xplot;
uexplot=eval(uexact);
if length(uexplot)==1 uexplot=ones(size(x))*uexplot; end
subplot(311),plot(xplot,uexplot,':')
figure(gcf)
title('exact solution')
x=xcopy;
end

%
%main loop
%
disp('entering main loop')
ok=-1;
%run=0;
xold=[xmin,xmax];
uold=[0,0];
while ok<0
%
%creating data vectors
%
xcopy=x;
h=diff(x);

x=xcopy(1:ni)+h/2;
onessizex=ones(size(x));
d=eval(diffusion);
if length(d)==1 d=onessizex*d; end
	if d < zeros(size(x))
         disp( 'error: d has to be positive'), return, end
%end	

c=eval(convection);
if length(c)==1 c=onessizex*c; end
r=eval(absorbtion);
if length(r)==1 r=onessizex*r; end
f=eval(force);
if length(f)==1 f=onessizex*f; end
x=x+0.01*h;
dplus=eval(diffusion);
x=x-0.01*h;
dminus=eval(diffusion);
Dd=(dplus-dminus)./(0.02*h);
clear onessizex; 
%
%plotting current mesh
%
x=xcopy;
hh=[h(1),h,h(ni)];
hh=(hh(2:ni+2)+hh(1:ni+1))/2;
subplot(313),plot(x,hh)
figure(gcf)
if run==0;
title('initial meshsize');
elseif run>0
title('new meshsize')
end
%pause(1);
%
%computing stiffness matrix (diagonal, left and right subdiagonals) and r.h.s.
%
[ad,al,ar]=stiff(ni,x,d,Dd,c,r,k,bctype);
b=loads(ni,x,d,Dd,c,r,f,g,bctype);
%
%solving tridiagonal system of equations
%
ndf=length(ad);
%ndf=ni-1+sum(bctype);
A=spdiags([ [al 0]' ad' [0 ar]' ],-1:1,ndf,ndf);
disp('starting to solve tridiagonal system')
u=(A\b')';
disp('finished solving tridiagonal system')
clear A;
%
%adding given boundary values
%
if bctype(1)==0 u=[g(1),u]; end
if bctype(2)==0 u=[u,g(2)]; end
%
length_of_u=length(u)
%
%plotting current solution
%
subplot(313),title('current meshsize')
if exactflag=='y'
subplot(311),plot(xplot,uexplot,':'),hold on,title('exact solution');
%pause(1)
end
%subplot(313),title('meshsize')
subplot(311),plot(x,u,'y-')
title('current solution')
%pause(2);
hold off
%
%testing if solution ok
%
[ok,globalerr,res]=test(ni,x,u,xold,uold,d,Dd,c,r,f,errnorm,tol,ad,al,ar,bctype);
%
%refining the mesh if error too large
%
if ok<0
disp('error too large')
error=globalerr
disp('refining the mesh')
%
if exactflag=='y';
exacterror=exacterr(ni,x,u,uexact,errnorm)
end
%
%pause(1)
%
%[ni,x]=corr(ni,x,res,errnorm,tol);
[ni,x]=corrnew(ni,x,res,errnorm,tol);
run=run+1;
%
%accepting current solution if error ok
%
elseif ok>0
disp('error ok');
error=globalerr
%
if exactflag=='y';
exacterror=exacterr(ni,x,u,uexact,errnorm)
if errnorm=='en'&exacterror<0.1*globalerr;
disp('I notice that the exact error is much smaller than the estimated one')
disp('Probably the (midpoint) quadrature rule has underestimated')
disp('the true error !')
end
if exacterror>globalerr;
disp('I notice that the exact error is larger than the estimated one')
disp('Please check that your exact solution is correctly specified')
disp('If so, please notify the author of the code')
end
end
%
%pause(1)
subplot(311),plot(x,u,'g-')
title('tested and accepted solution')
%title('solution')
%pause(1);
end %test loop
end %main loop

------------new file--------------

function [ni,x]=corr(ni,x,res,errnorm,tol)
disp('corr')
h=diff(x);
%
%computing new mesh size
%
if errnorm=='en';
hnew=(tol^2./(res+eps).^2/ni).^(1/3);
%
elseif errnorm=='l2'
hnew=(tol^2./(res+eps).^2/ni).^(1/5);
%
elseif errnorm=='ma'
hnew=sqrt(tol./(res+eps));
end
%
%computing new mesh points
%
%sons=max(ones(size(h)),ceil(h./hnew));
sons=min(8,max(ones(size(h)),ceil(h./hnew))); %limit the number of sons
%totalsum=sum(sons);
%maxallowed = 50000;
hnew=h./sons;
cumsum(sons);
xnew=x(1);
for j=1:ni;
add=[x(j)+hnew(j):hnew(j):x(j+1)];
xnew=[xnew,add];
end
x=xnew;
ni=length(x)-1;
for i=1:3
x=[x(1), 0.5*(x(1:ni-1)+x(3:ni+1)), x(ni+1)]; %smoothing of new grid
end

----------------new file---------------

function [ni,x]=corrnew(ni,x,res,errnorm,tol)
disp('corrnew')
h=diff(x);
%
%computing new mesh size
%
if errnorm=='en';
hnew=(tol^2./(res+eps).^2/ni).^(1/3);
%
elseif errnorm=='l2'
hnew=(tol^2./(res+eps).^2/ni).^(1/5);
%
elseif errnorm=='ma'
hnew=sqrt(tol./(res+eps));
end
%
%computing new mesh points
%
sons=min(12,max(ones(size(h)),ceil(h./hnew))); %limit the number of sons
cumulsum=cumsum([ 1 sons(1:length(sons)-1)]);
%cumulsum=cumulsum(1:length(cumulsum)-1);
maxson=max(sons);
hnew=h./sons;
xnew(cumulsum)=x(1:ni);
for i=1:maxson-1
xnew(cumulsum+min(i,sons)) = x(1:ni) + min(i,sons).*hnew;
end
if xnew(length(xnew)) ~= x(ni+1)
  xnew=[xnew x(ni+1)];
end
x=xnew;
ni=length(x)-1;
for i=1:6
x=[x(1), 0.5*(x(1:ni-1)+x(3:ni+1)), x(ni+1)]; %smoothing of new grid
end


-------------new file----------

function b=loads(ni,x,d,Dd,c,r,f,g,bctype);
disp('load')
%h=x(2:ni+1)-x(1:ni);
h=diff(x);
%
%setting (local) constant in streamline diffusion test quantity: v+chv'
%
csd=sign(c)/2.*( abs(c).*h > d);
%
j=1:ni-1;
b(j)=f(j).*(0.5+csd(j)).*h(j)+f(j+1).*(0.5-csd(j+1)).*h(j+1);
%
if bctype(1)==0
b(1)=b(1)-(-d(1)/h(1)-c(1)/2+r(1)*h(1)/6)*g(1);
b(1)=b(1)-csd(1)*(Dd(1)-c(1)+r(1)*h(1)/2)*g(1);
else
b=[f(1)*(1/2-csd(1))*h(1)+g(1),b];
end
if bctype(2)==0
ndf=ni-1+bctype(1);
b(ndf)=b(ndf)-(-d(ni)/h(ni)+c(ni)/2+r(ni)*h(ni)/6)*g(2);
b(ndf)=b(ndf)-csd(ni)*(Dd(ni)-c(ni)-r(ni)*h(ni)/2)*g(2);
else
b=[b,f(ni)*(1/2+csd(ni))*h(ni)+g(2)];
end


----------------new file-----------------

function [ad,al,ar]=stiff(ni,x,d,Dd,c,r,k,bctype)
disp('stiff')
h=diff(x);
%
%setting constant in streamline diffusion test quantity: v+chv'
%
csd=sign(c)/2.*(abs(c).*h>d);
if (csd > 0);
disp('using streamline diffusion modification');
end
%

i=2:ni-1;
j=2:ni-2;
ad(i-1)=d(i)./h(i)-c(i)/2+r(i).*h(i)/3+csd(i).*(-Dd(i)+c(i)-r(i).*h(i)/2);
ad(j)  = ad(j)+d(j)./h(j)+c(j)/2+r(j).*h(j)/3+csd(j).*(-Dd(j)+c(j)+r(j).*h(j)/2);
j=ni-1;
ad(j)=d(j)./h(j)+c(j)/2+r(j).*h(j)/3+csd(j).*(-Dd(j)+c(j)+r(j).*h(j)/2);
ad(1)=ad(1)+d(1)/h(1)+c(1)/2+r(1)*h(1)/3+csd(1)*(-Dd(1)+c(1)+r(1)*h(1)/2);
al(i-1)=-d(i)./h(i)-c(i)/2+r(i).*h(i)/6+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+csd(i).*(Dd(i)-c(i)-r(i).*h(i)/2);

ad(ni-1)=ad(ni-1)+d(ni)/h(ni)-c(ni)/2+r(ni)*h(ni)/3+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)];
ar=[ar,-d(ni)/h(ni)+c(ni)/2+r(ni)*h(ni)/6+csd(ni)*(Dd(ni)-c(ni)-r(ni)*h(ni)/2)];
end

----------------new file-------------

function [ni,x]=pred(xmin,xmax,f,errnorm,tol)
disp('pred')
l=xmax-xmin;
l2=l*l;
ni=max(5,round(sqrt(0.2*l/tol)));
if errnorm=='en';
ni=ni*ni;
end
ni=min(20,ni);  % limit maximal size
dx=l/ni;
x=[xmin:dx:xmax];

-----------------new file----------------

function [ok,globalerr,res]=test(ni,x,u,xold,uold,d,Dd,c,r,f,errnorm,tol,ad,al,ar,bctype)
disp('test')
h=diff(x);
%
i=1:ni;
res=f+(Dd-c).*(u(i+1)-u(i))./h-r.*(u(i+1)+u(i))/2;
%
%plotres=[res(1),res,res(ni)];
%plotres=(plotres(1:ni+1)+plotres(2:ni+2))/2;
%subplot(312),plot(x,plotres)
%title('residual')
%pause(3)
%
if errnorm=='en';
%
res=abs(res)./sqrt(d);
%
cint=1/sqrt(10); %interpolation error constant in ||z-zi||<c||z'|| (where z=e)
%
%computing global energy norm error
%
res=cint*res;
hres=h.*res;
globalerr=sqrt(sum(hres.*hres.*h));
cstab=1
%
elseif errnorm=='l2';
%
xcopy=x;
niold=size(xold,2)-1;
ndf=ni-1+sum(bctype);
%
cint=1/sqrt(30)/2;    %interpolation error constant
%
%setting up dual problem
%


x=xcopy(1)+h(1)/2;
wl=(xold(2)-x)/(xold(2)-xold(1));
wr=(x-xold(1))/(xold(2)-xold(1));
%%%ee=(u(2)+u(1))/2-(wl*uold(1)+wr*uold(2));
ee=rand;
eefirst=ee;
l2e=ee*ee*h(1);
e(1)=ee*h(1)/2;
index=1;

for i=2:ni-1;
x=xcopy(i)+h(i)/2;

while xold(index)<x;
index=index+1;
end

wl=(xold(index)-x)/(xold(index)-xold(index-1));
wr=(x-xold(index-1))/(xold(index)-xold(index-1));
%%%ee=(u(i+1)+u(i))/2-(wl*uold(index-1)+wr*uold(index));
ee=rand;
l2e=l2e+ee*ee*h(i);
e(i-1)=e(i-1)+ee*h(i)/2;
e(i)=ee*h(i)/2;
end

x=xcopy(ni)+h(ni)/2;
wl=(xold(niold+1)-x)/(xold(niold+1)-xold(niold));
wr=(x-xold(niold))/(xold(niold+1)-xold(niold));
%%%ee=(u(ni)+u(ni+1))/2-(wl*uold(niold)+wr*uold(niold+1));
ee=rand;
eelast=ee;
l2e=l2e+ee*ee*h(ni);
l2e=sqrt(l2e);
e(ni-1)=e(ni-1)+ee*h(ni)/2;
%
if bctype(1)==1;
e=[eefirst*h(1)/2,e];
end
if bctype(2)==1;
e=[e,eelast*h(ni)/2];
end
%
%solving dual problem
%
%
%adcopy=ad;
%for i=2:ndf;
%ad(i)=ad(i)-al(i-1)*ar(i-1)/ad(i-1);
%e(i)=e(i)-e(i-1)*ar(i-1)/ad(i-1);
%end
%z(ndf)=e(ndf)/ad(ndf);
%for i=ndf-1:-1:1;
%e(i)=e(i)-e(i+1)*al(i)/ad(i+1);
%z(i)=e(i)/ad(i);
%end
%ad=adcopy;
%
ndf=length(ad);
%ndf=ni-1+sum(bctype);
A=spdiags([ [ar 0]' ad' [0 al]' ],-1:1,ndf,ndf);  % note ar and al change places
z=(A\e')';
clear A;
%
if bctype(1)==0
z=[0,z];    %given boundary value
end
if bctype(2)==0
z=[z,0];    %given boundary value
end

%
%computing D2z and cstab
%
Dz=(z(2:ni+1)-z(1:ni))./h;
DDz=abs(Dz(2:ni)-Dz(1:ni-1));
DDz=[DDz(1),DDz,DDz(ni-1)];
DDz=(DDz(2:ni+1)+DDz(1:ni))/2;
DDz=DDz./h;
cstab=sqrt(sum(DDz.*DDz.*h))/l2e
%
res=cstab*cint*abs(res);
hhres=h.*h.*res;
globalerr=sqrt(sum(hhres.*hhres.*h));
%
x=xcopy;
%
elseif errnorm=='ma';
%
ndf=ni-1+sum(bctype);
%
cint=1/8; %interpolation constant
%
%loop over test points to determine cstab
%
A=spdiags([ [ar 0]' ad' [0 al]' ],-1:1,length(ad),length(ad)); % ar and al change places
%
cstab=0;
cstabdiff=0;
globalerr=0;
nitp=5; %number of internal points
index=[1:max(1.01,ndf/(nitp+1)):ndf];
index=round(index);
index=[index,ndf];
si=size(index,2);
%
delta=sparse(zeros(ndf,si));
%
%setting up and solving dual problem(s)
%
for j=1:si  delta(index(j),j)=1;  end
%
zz=(A\delta)';
%
if bctype(1)==0
zz=[zeros(si,1),zz];    %given boundary value
end
if bctype(2)==0
zz=[zz,zeros(si,1)];    %given boundary value
end
%
for j=1:si
z=full(zz(j,:));  % using one of the right hand sides
%
%computing D2z and cstab
%
Dz=(z(2:ni+1)-z(1:ni))./h;
DDz=abs(Dz(2:ni)-Dz(1:ni-1));
%DDz=(Dz(2:ni)-Dz(1:ni-1));
DDz=[DDz(1),DDz,DDz(ni-1)];
DDz=(DDz(2:ni+1)+DDz(1:ni))/2;
DDz=DDz./h;
cstab=max(cstab,sum(abs(DDz).*h));
%
Duu=(u(2:ni+1)-u(1:ni))./h;
DDuu=(Duu(2:ni)-Duu(1:ni-1));
DDuu=[DDuu(1),DDuu,DDuu(ni-1)];
DDuu=(DDuu(2:ni+1)+DDuu(1:ni))/2;
DDuu=DDuu./h;
csd=sign(c)/2.*(abs(c).*h>d);
differr=abs(sum(csd.*h.^2.*(d.*DDuu(1:size(h,2))+Dd.*Duu).*Dz(1:size(h,2))));
cstabdiff=max(cstabdiff,abs(sum(abs(csd).*h.*abs(Dz(1:size(h,2))))));
if globalerr < differr+abs(cint*sum(DDz.*(res).*h.^3))
	globalerr=max(globalerr,differr+abs(cint*sum(DDz.*(res).*h.^3)));
	streamline_error=differr;
end
%
end % for j=1:si
%
cstab
streamline_error
%globalerr
res=cstab*cint*abs(res)+cstabdiff*(d.*abs(DDuu)+Dd.*Duu)./h;
%cstab
%res=cstab*cint*abs(res);
%globalerr=max((h.^2).*res);
%
end %of errnorm branching 
%
%subplot(312),title(['residual, stability factor =',num2str(cstab)])
%
%setting ok flag
%
if globalerr<tol
ok=1;
else
ok=-1;
end

---------------------new file---------------------


p=1;
qrule=1;
%set(1,'PaperUnits','centimeters')
%set(1,'PaperPosition',[5,5,10,12])
%set(1,'PaperType','a4letter')


---------------------new file---------------------

function exacterror=exacterr(ni,x,u,uexact,errnorm)
global qrule
%gauss=1/sqrt(3);
%gauss1=1/2-gauss/2;
%gauss2=1/2+gauss/2;
ucomp=(u(2:ni+1)+u(1:ni))/2;
%ucomp(1:2:2*ni-1)=(gauss1*u(2:ni+1)+(1-gauss1)*u(1:ni));
%ucomp(2:2:2*ni)=(gauss2*u(2:ni+1)+(1-gauss2)*u(1:ni));
%h=x(2:ni+1)-x(1:ni);
h=diff(x);
Ducomp=(u(2:ni+1)-u(1:ni))./h;
%Ducomp(1:2:2*ni-1)=(u(2:ni+1)-u(1:ni))./h;
%Ducomp(2:2:2*ni)=Ducomp(1:2:2*ni-1);
xcopy=x;
%

%for i=1:ni;
%x=xcopy(i)+h(i)/2;
%%x=xcopy(i)+gauss1*h(i);
%uex(i)=eval(uexact);
%%uex(2*i-1)=eval(uexact);
%%x=xcopy(i)+qrule*h(i)/2;
%Duex(i)=eval(uexact);
%%Duex(2*i-1)=eval(uexact);
%x=x+0.01*h(i);
%Duex(i)=(eval(uexact)-Duex(i))/(0.01*h(i));
%%Duex(2*i-1)=(eval(uexact)-Duex(2*i-1))/(0.01*h(i));
%%
%%x=xcopy(i)+gauss2*h(i);
%%uex(2*i)=eval(uexact);
%%Duex(2*i)=eval(uexact);
%%x=x+0.01*h(i);
%%Duex(2*i)=(eval(uexact)-Duex(2*i))/(0.01*h(i));
%end

x=xcopy(1:ni)+h/2;
uex=eval(uexact);
if length(uex)==1 uex=ones(size(x))*uex; end
x=x+0.01*h;
uexplus=eval(uexact);
if length(uexplus)==1 uexplus=ones(size(x))*uexplus; end
Duex=(uexplus-uex)./(0.01*h);

%
%hh(1:2:2*ni-1)=h;
%hh(2:2:2*ni)=h;
%abs(Ducomp-Duex)
if errnorm=='en';
exacterror=sqrt(sum((Ducomp-Duex).^2.*h));
elseif errnorm=='l2';
exacterror=sqrt(sum((ucomp-uex).^2.*h));
elseif errnorm=='ma';
exacterror=max(abs(uex-ucomp));
end

--------------end of file-----------------