%w1,x1,w2,x2
f=inline('[1-w1-w2;1/2-w1*x1-w2*x2;1/3-w1*x1^2-w2*x2^2;1/4-w1*x1^3-w2*x2^3]',
'x1','w1','x2','w2');
df=inline('[0,-1,0,-1;-w1,-x1,-w2,-x2;-2*w1*x1,-x1^2,-2*w2*x2,-x2^2;-3*w1*x1^2,-x1^3,-3*w2*x2^2,-x2^3]',
'x1','w1','x2','w2');

u=[0,1,1,1];
eps=1e-8;
for i=1:100
	x1=u(1);w1=u(2);x2=u(3);w2=u(4);
	u=[x1;w1;x2;w2]-df(x1,w1,x2,w2)\f(x1,w1,x2,w2);
	if abs(u-[x1;w1;x2;w2])<eps,break,end;
end;

x1=u(1);w1=u(2);x2=u(3);w2=u(4);
g=inline('exp(-x.^2)','x');
n=10;
a=0; b=1;
h=(b-a)/n;

I=h*sum([w1*g(h*((0:(n-1))+x1)),w2*g(h*((0:(n-1))+x2))]);
printf('I=%0.10f\n',I);

f=inline('x.^2+y.^2','x','y');
n=10; h=1/n; eps=1e-8;
[x,y]=meshgrid(0:h:1,0:h:1); U=f(x,y);
U(2:end-1,2:end-1)=0; mesh(U);
%print -dtex slika13-1.eps
for k=1:1000
  V=U;
  for i=2:n,
   for j=2:n,
     U(i,j)=(U(i+1,j)+U(i-1,j)+U(i,j-1)+U(i,j+1))/4;
   end;
  end;
  if abs(U-V)<eps, break, end;  
end;
figure; mesh(U);
%print -dtex slika13-2.eps
