  
diffusion[n_,type0_,lim0_,type1_,lim1_,coef_,vecS_,vecU_,temps_,dt_,tf_]:=Module[
    {t,u,dx,matA,matB,vecC,vecD,a,nd,alpha,j1,j2,k,j},
    t = temps;
    u = vecU;
    dx=1/(n-1);
    matA = SparseArray[{},{n,n}];
    matB = SparseArray[{},{n,n}];
    vecC = vecS*dt;
    a = 2*dt/dx^2;
    nd = Length[coef];
    vecD = Array[0,{n,1}];
    j1=1;
    For[k=1,k<=nd,k++,
        j2 = Floor[coef[[k]][[1]]*n];
        d = coef[[k]][[2]];
        For[j=j1,j<=j2,j++,
            vecD[[j]] = d
        ];
        j1 = j2;
    ]; 
    For[j=2,j<=n-1,j++,
        alpha = a*vecD[[j]];
        matA[[j,j-1]] = -alpha/2; matB[[j,j-1]] = alpha/2;
        matA[[j,j]] = 1+alpha; matB[[j,j]] = 1-alpha;
        matA[[j,j+1]] = -alpha/2; matB[[j,j+1]] = alpha/2; 
    ];
    If[type0=="dirichlet",
        matA[[1,1]] = 1; matA[[1,2]] = 0; vecC[[1]] = lim0;
    ];
    If[type0=="neumann",
        matA[[1,1]] = -1; matA[[1,2]] = 1; vecC[[1]] = dx*lim0;
    ];
    If[type1=="dirichlet",
        matA[[n,n-1]] = 0; matA[[n,n]] = 1; vecC[[n]] = lim1;
    ];
    If[type1=="neumann",
        matA[[n,n-1]] = 1; matA[[n,n]] = -1; vecC[[n]] = dx*lim1;
    ];
    For[k=1,k<=nd-1,k++,
        (*frontiere entre deux coefficients de diffusion differents*)
        j = Floor[coef[[k]][[1]]*n]-1; 
        matA[[j,j-1]] = -vecD[[j]]; 
        matA[[j,j]] = vecD[[j]]+vecD[[j+1]]; 
        matA[[j,j+1]] = -vecD[[j+1]];
        matB[[j,j-1]] = 0; matB[[j,j]] = 0; matB[[j,j+1]] = 0; vecC[[j]] = 0;
    ];
    While[t<tf,
        u = LinearSolve[matA,matB.u+vecC];
        t = t+dt;
    ];
    Return[{u,t}];
]
            