function dE=denergy1(vdoft)
global A
global L B pp nel
global npe nodg ndof conecg conecdof coordg
vdof=vdoft';
xh=-1:0.1:1;#integration points
wh=[1/20 1/10*ones(1,19) 1/20];#integr. weights
nhat=size(xh,2);
##[g dg d2g] = master(1,xh);
[g dg d2g]=masterlag(1,xh);
[sf dsf ddsf]=masterlag(pp,xh);
dE=zeros(ndof,1);
for k=1:nel
 XX(:,1)=coordg(:,conecg(k,1));
 XX(:,2)=coordg(:,conecg(k,2));
 x=zeros(1,nhat);
 dx=zeros(1,nhat);
 for i=1:2
  x(1,:)=x(1,:)+XX(1,i)*g(i,:);
  dx=dx+XX(1,i)*dg(i,:);
 end
 uh=vdof(conecdof(k,:))*sf;
 duh=vdof(conecdof(k,:))*dsf./dx;
####data
 beta=ones(1,nhat)+4*(1+sign(x-2*L/3)).*ones(1,nhat);
 gamma=2*ones(1,nhat);
 ff=x.*x.*ones(1,nhat);
####
 for i=1:pp+1
  dE(conecdof(k,i))=dE(conecdof(k,i))+...
     (dsf(i,:)./dx)*(duh.*beta.*wh.*dx)'+...
     sf(i,:)*(uh.*gamma.*wh.*dx)'-...
     sf(i,:)*(ff.*wh.*dx)';
 end
end
dE(ndof)=dE(ndof)-B;
dE(1)=vdof(1)-A;
end