function [z]=f(x);
    z=exp(-(x-1)^2/2)/sqrt(2*%pi);
endfunction

function [z]=g(x)
    z=1/(1+x^2);
endfunction

function [z]=h(x)
    z=f(x)*g(x);
endfunction

function [z]=transi(x);
    z=0.9*x+grand(1,1,'nor',0,1);
endfunction

function [z]=M(x,y)
    z=exp(-(y-0.9*x)^2)/sqrt(2*%pi);
endfunction

function [z]=ratio(x,y)
    z=(M(y,x)*f(y)*g(y))/(M(x,y)*f(x)*g(x));
endfunction

N=10000;
x=0;
valeurs=[x];
for i=1:N
    //disp("coucou");
    y=transi(x);
    u=grand(1,1,'unf',0,1);
    //disp(ratio(x,y));
    //disp(x);
    if (u<ratio(x,y)) then
        x=y;
    end
    valeurs=[valeurs, x];
end
clf()
histplot(ceil(N/100),valeurs,style=2);
u=[-3:0.1:3];
z=feval(u,h);
C=integrate('f(x)*g(x)','x',-3,5);
plot(u,z/C,'red');
