function traj = heun1d(mu,sigma,x0,h,n,ntraj)
//      Retourne une matrice a n+1 lignes et nbtraj colonnes.
//      Les colonnes contiennent des simulations independantes 
//      de trajectoires du processus de diffusion de derive mu 
//      et diffusion sigma (fonctions de deux variables).
//      Les trajectoires partent de x0 et sont simulees par la methode
//      d'Euler de pas h sur n pas.

sh = sqrt(h);
pas_brownien = grand(n,ntraj,"nor",0,sh);
deff("y=mu_c(t,x)","y=mu(t,x)-sigma(t,x).^2/2");
traj = zeros(n+1,ntraj);
t = 0; X = x0*ones(1,ntraj);
traj(1,:) = X;
for i=1:n,
    Xeuler = X + mu_c(t,X)*h + sigma(t,X).*pas_brownien(i,:);
    X = X + (mu_c(t,X)+mu_c(t+h,Xeuler))*h/2 +..
            (sigma(t,X)+sigma(t+h,X)).*pas_brownien(i,:); 
    traj(i+1,:) = X;
    t = t+h;
end;
xbasc();
tmin=0; tmax=n*h;
xmin = min(traj); xmax=max(traj)
plotframe([tmin,xmin,tmax,xmax],[10,1,10,1],[%f,%f],["Trajectoires de diffusion","",""]);
t=[0:h:n*h]'*ones(1,nbtraj);
plot2d(t,traj,ones(1,ntraj),"000");



