u = 1 //velocidade de propagação da onda
p = 1 //periodo
N = 100//valor
te=0.8;
dt = 0.008 //variação do tempo
L = 1;
dx = L/N; //variação do espaço
C = u*(dt/dx) //número de Courant
passodetempo=tempomax/dt
count=0;
x = 0:dx:dx
tempomax=0.5
c0=0
condcontorno=1 //concentração
Q0=zeros(1,N+5)
//solanalitica
solanalitica=Q0
i=1
for j=0:dt*2:te
    solanalitica(i)=condcontorno
    i=i+1
end
disp (solanalitica)
function superbee = superbe(theta)
        superbee = max(0,max (min(1,2*theta),min(2,theta)));
endfunction

function mc = mcc(theta)
        mc= max(0,min(min((1+theta)/2,2),2*theta ) );
endfunction

function vanleer=vanlee(theta)
         vanleer = (theta+abs(theta))/(1+abs(theta));
endfunction

function thetanegativo= tneg(Q, i)

    if((Q(i) - Q(i-1)) == 0) then
        thetanegativo= 0;
    else
        thetanegativo= (Q(i-1) - Q(i-2))/(Q(i) - Q(i-1));
    end

endfunction

function thetapositivo= tpos(Q, i)

    if ((Q(i+1)-Q(i)) == 0) then

        thetapositivo= 0;

    else

        thetapositivo= (Q(i)-Q(i-1))/(Q(i+1) - Q(i));

    end

endfunction

//botando c0 no vetor de zeros e a cond de contorno
//atribuindo o a condição inicial no vetos a partir da segunda posicao
atrQmaisumUpwind=zeros(1,N+5)
atrQUpwind=zeros(1,N+5)
atrQmaisumVan=zeros(1,N+5)
atrQVan=zeros(1,N+5)
atrQmaisumMC=zeros(1,N+5)
atrQMC=zeros(1,N+5)
atrQmaisumSuperB=zeros(1,N+5)
atrQSuperB=zeros(1,N+5)

for i = 1 : 1
    atrQmaisumUpwind(i) = condcontorno;
    atrQUpwind(i)=condcontorno;
    atrQmaisumVan(i) = condcontorno;
    atrQVan(i) = condcontorno;
    atrQmaisumMC(i)=condcontorno;
    atrQMC(i)=condcontorno;
    atrQmaisumSuperB(i)=condcontorno;
    atrQSuperB(i)=condcontorno;
    Q0(i)=condcontorno;
end
for i=2:N+5
    atrQmaisumUpwind(i)=c0;
    atrQUpwind(i)=c0;
    atrQmaisumVan(i) = c0;
    atrQVan(i) = c0;
    atrQmaisumMC(i)=c0;
    atrQMC(i)=c0;
    atrQmaisumSuperB(i)=c0;
    atrQSuperB(i)=c0;
    Q0(i)=c0;
end
//disp(atrQSuperB)

//utilizando a expressão fornecida
for n = 1 : passodetempo

    i=2;
//utilizando o metodo upwind pra atribuir o valor e começar as contas
    atrQmaisumUpwind(i) = atrQUpwind(i) - C*(atrQUpwind(i)-atrQUpwind(i-1));
    atrQmaisumSuperB(i) = atrQSuperB(i) - C*(atrQSuperB(i)-atrQSuperB(i-1));
    atrQmaisumVan(i) = atrQVan(i) - C*(atrQVan(i)-atrQVan(i-1));
    atrQmaisumMC(i) = atrQMC(i) - C*(atrQMC(i)-atrQMC(i-1));

        for i = 3 : N
            //upwind
            atrQmaisumUpwind(i) = atrQUpwind(i) - C*(atrQUpwind(i)-atrQUpwind(i-1));
            //vanlee
            D=(atrQVan(i) - C*(atrQVan(i)-atrQVan(i-1)))
            F=((1/2)*C*(1 - C))*(vanlee(tpos(atrQVan,i))*(atrQVan(i+1)-atrQVan(i))-vanlee(tneg(atrQVan, i))*(atrQVan(i)-atrQVan(i-1)))
            atrQmaisumVan(i) =  (D)-(F)


            //superbee
            A= (atrQSuperB(i) - C*(atrQSuperB(i)-atrQSuperB(i-1)))
            B=((1/2)*C*(1 - C))*(superbe(tpos(atrQSuperB,i))*(atrQSuperB(i+1)-atrQSuperB(i))-superbe(tneg(atrQSuperB, i))*(atrQSuperB(i)-atrQSuperB(i-1)))
            atrQmaisumSuperB(i) =(A)-(B)

            //mc
            G=atrQMC(i) - C*(atrQMC(i)-atrQMC(i-1))
            H=((1/2)*C*(1-C))*(mcc(tpos(atrQMC, i))*(atrQMC(i+1)-atrQMC(i)) - mcc(tneg(atrQMC, i))*(atrQMC(i)-atrQMC(i-1)))
            atrQmaisumMC(i) = (G)-(H)
            count = count + 1;
        end

    x = x + dt;
    atrQSuperB = atrQmaisumSuperB;
    atrQMC = atrQmaisumMC;
    atrQUpwind = atrQmaisumUpwind;
    atrQVan = atrQmaisumVan;
end

figure(1)
titulo='Solução para Superbee para t=0.8'
xtitle(titulo)
plot(atrQSuperB,'y-','thickness', 2)
plot(solanalitica,'b-','linest','--')
a=get("current_axes")
a.data_bounds=[-10,-0.5;110,1.5];

figure(2)
titulo='Solução Van Leer para t=0.8'
xtitle(titulo)
plot(atrQVan,'y-','thickness', 2)
plot(solanalitica,'b-','linest','--')
a=get("current_axes")
a.data_bounds=[-10,-0.5;110,1.5];

figure(3)
titulo='Solução Upwind para t=0.8'
xtitle(titulo)
plot(atrQUpwind,'y-','thickness', 2)
plot(solanalitica,'b-','linest','--')
a=get("current_axes")
a.data_bounds=[-10,-0.5;110,1.5];

figure (4)
titulo='Solução MC para t=0.8'
xtitle(titulo)
plot(atrQMC,'y-','thickness', 2)
plot(solanalitica,'b-','linest','--')
a=get("current_axes")
a.data_bounds=[-10,-0.5;110,1.5];