
function [U,t]=diffusion(N,type0,lim0,type1,lim1,coef,S,U,t,dt,tf)
    dx=1/(N-1)
    A=spzeros(N,N)
    B=spzeros(N,N)
    C=S*dt
    a=2*dt/dx^2
    nd=size(coef)
    nd=nd(1)
    D = zeros(N) // coefficients de diffusion
    j1=1
    for k=1:nd,
        j2 = int(coef(k,1)*N)
        d = coef(k,2)
        for j=j1:j2,
            D(j)=d
        end
        j1 = j2
    end
    for j=2:N-1, // schema de Cranck-Nicholson
        A(j,j-1)=-a*D(j)/2; B(j,j-1)=a*D(j)/2;
        A(j,j)=1+a*D(j); B(j,j)=1-a*D(j);
        A(j,j+1)=-a*D(j)/2; B(j,j+1)=a*D(j)/2;
    end;
    if type0=="dirichlet" then
        A(1,1)=1;
        A(1,2)=0;
        C(1)=lim0;
    end
    if type0=="neumann" then
        A(1,1)=-1;
        A(1,2)=1;
        C(1)=dx*lim0;
    end
    if type1=="dirichlet" then
        A(N,N-1)=0;
        A(N,N)=1;
        C(N)=lim1;
    end
    if type1=="neumann" then
        A(N,N)=1;
        A(N,N-1)=-1;
        C(N)=dx*lim1;
    end
    for k=1:nd-1,
        j = int(coef(k,1)*N)-1 // frontiere entre deux coefficients de diffusion differents
        A(j,j-1)=-D(j); A(j,j)=D(j)+D(j+1); A(j,j+1)=-D(j+1);
        B(j,j-1)=0; B(j,j)=0; B(j,j+1)=0; C(j)=0;
    end
    while t<tf // boucle de calcul
        U = lusolve(A,B*U+C);
        t = t+dt
    end
endfunction
             