xbasc();
n=10000;
m1=0;
m2=3;
function [y]=f(x)
  y=exp(-(x-m1)^2/2)/sqrt(2*%pi);
endfunction
function [y]=g(x)
  y=exp(-(x-m2)^2/2)/sqrt(2*%pi);
endfunction
a=zeros(n,1);
for i=1:n
  b=0;
  while (b==0)
    u=rand();
    v=rand();
    x=grand(1,1,'nor',m1,1);
    y=grand(1,1,'nor',m2,1);
    if (u*f(x)<g(x)) then
      z=y;
      b=1;
    else 
      if (v*g(y)<f(y)) then
        z=x;
        b=1;
      end,
    end,
    a(i,1)=z;
  end,  
end
histplot(ceil(sqrt(n)),a,style=2); 
//on dessine aussi le graphique de h
t=grand(n,1,'unf',-5,5);
h=max(f(t),g(t));
s=10*sum(h)/n;
t=-5:0.1:5;
h=max(f(t),g(t))./s;
plot(t,h,'red');
