function [x,T]=tdma(u,N,S)
L=1;
rho=1;
r=0.1;
dx=L/N;
phia=4;
phib=1;

F=rho*u;
D=r/dx;
Pe=F/D;

b(1)=1;
c(1)=0;
d(1)=phia;

switch S
    case 'central'
        %boundary near the base,x=dx/2
        aw=0;
        ae=D-F/2;
        su=(2*D+F)*phia;
        sp=-(2*D+F);
        ap=aw+ae-sp;
        
        a(1)=-aw;
        b(2)=ap;
        c(2)=-ae;
        d(2)=su;
        
        %inner points
        aw=D+F/2;
        ae=D-F/2;
        su=0;
        sp=0;
        ap=aw+ae-sp;
        
        for i=3:N
            a(i-1)=-aw;
            b(i)=ap;
            c(i)=-ae;
            d(i)=su;
        end
        
        %end boundary, at x=L
        aw=D+F/2;
        ae=0;
        su=-(2*D-F);
        sp=(2*D+F)*phib;
        ap=aw+ae-sp;
        
        a(N)=-aw;
        c(N+1)=-ae;
        b(N+1)=ap;
        d(N+1)=su;
        
    case 'upwind'
        %boundary near the base,x=dx/2
        aw=0;
        ae=D;
        su=(2*D+F)*phia;
        sp=-(2*D+F);
        ap=aw+ae-sp;
        
        a(1)=-aw;
        b(2)=ap;
        c(2)=-ae;
        d(2)=su;
        
        %inner points
        aw=D+max([F,0]);
        ae=D+max([0,-F]);
        su=0;
        sp=0;
        ap=aw+ae-sp;
        
        for i=3:N
            a(i-1)=-aw;
            b(i)=ap;
            c(i)=-ae;
            d(i)=su;
        end
        
        %end boundary, at x=L
        aw=D+F;
        ae=0;
        su=2*D*phib;
        sp=-2*D;
        ap=aw+ae-sp;
        
        a(N)=-aw;
        c(N+1)=-ae;
        b(N+1)=ap;
        d(N+1)=su;
        
    case 'hybrid'
        %boundary near the base,x=dx/2
        aw=0;
        ae=0;
        su=(2*D+F)*phia;
        sp=-(2*D+F);
        ap=aw+ae-sp;
        
        a(1)=-aw;
        b(2)=ap;
        c(2)=-ae;
        d(2)=su;
        
        %inner points
        aw=max([F,D+F/2,0]);
        ae=max([-F,D-F/2,0]);
        su=0;
        sp=0;
        ap=aw+ae-sp;
        
        for i=3:N
            a(i-1)=-aw;
            b(i)=ap;
            c(i)=-ae;
            d(i)=su;
        end
        
        %end boundary, at x=L
        aw=F;
        ae=0;
        su=2*D*phib;
        sp=-2*D;
        ap=aw+ae-sp;
        
        a(N)=-aw;
        c(N+1)=-ae;
        b(N+1)=ap;
        d(N+1)=su;
        
    case 'powlaw'
        %boundary near the base,x=dx/2
        aw=0;
        ae=0;
        su=(2*D+F)*phia;
        sp=-(2*D+F);
        ap=aw+ae-sp;
        
        a(1)=-aw;
        b(2)=ap;
        c(2)=-ae;
        d(2)=su;
        
        %inner points
        aw=D*max([0,(1-0.1*abs(Pe))^5])+max([F,0]);
        ae=D*max([0,(1-0.1*abs(Pe))^5])+max([-F,0]);
        su=0;
        sp=0;
        ap=aw+ae-sp;
        
        for i=3:N
            a(i-1)=-aw;
            b(i)=ap;
            c(i)=-ae;
            d(i)=su;
        end
        
        %end boundary, at x=L
        aw=F;
        ae=0;
        su=2*D*phib;
        sp=-2*D;
        ap=aw+ae-sp;
        
        a(N)=-aw;
        c(N+1)=-ae;
        b(N+1)=ap;
        d(N+1)=su;
end
a(N+1)=0;
b(N+2)=1;
d(N+2)=phib;
%start tdma
bd(1)=b(1);
dd(1)=d(1);

for i=2:N+2
    bd(i)=b(i)-(a(i-1)*c(i-1))/bd(i-1);
    dd(i)=d(i)-(a(i-1)*dd(i-1)/bd(i-1));
end
T(N+2)=dd(N+2)/bd(N+2);
for i=N+1:-1:1
    T(i)=(dd(i)-c(i)*T(i+1))/bd(i);
end
%end tdma
%matching x axis
x(1)=0;
x(2)=dx/2;
for i=3:N+1
    x(i)=x(i-1)+dx;
end
x(N+2)=L;