//(Scilab)
//The following code allows to solve the 1-dimentionnal heat equation on the half space \R_+ provided with Robin boundary conditions (Crank-Nicolson scheme)

//****: changeable parameters

clear()

//ROBIN COEFFICIENT:
//alpha=0.99;   //Dirichlet (do not take alpha=1 because one divides by 1-alpha somewhere bellow...)
//alpha=0;      //Neumann
alpha=0.5;      //Non-degenerated Robin cases (0<alpha<1)****

L=10;  //length of domain****
d=10;        //diffusion coefficient****
N=500;      //nb of descretisation points****
Tmax=0.1;    //final time****
u0max=1000;   //maximum of initial datum****

h=L/(N-1);      //computation of space step
cfl=0.8;
dt=cfl*(h*h)/(2*d);      //computation of time step

XX=0:h:L;   //space mesh
TT=0:dt:Tmax;   //time mesh
MM=zeros(TT);

function y=arr(x,n)//rounding function, eg.: arr(2.564,2)=2.56
    y=round(x*10^(n))*10^(-n);
endfunction

A=alpha/(1-alpha);  //alpha u-(1-alpha)u_x=0  <=> A u-u_x=0 (beware alpha is not 1)

lambda=(d*dt)/(h*h);    //used in matrix bellow

function y=u0(x)        //initial datum: indicator function of the ball B(1.5,0.5)
    if 1<x & x<2 then
    y=u0max;
    else
    y=0
    end
endfunction

function x=find_max_pos(A)
    x=XX(1);
    for i=2:length(A)
        if A(i)>A(i-1) then
        x=XX(i);
        end
    end
endfunction

U0=zeros(XX);
for i=1:length(XX)
    U0(i)=u0(XX(i));
end

Up=U0;  //previous U
U=zeros(XX);  //next U

//MATRIX A
matrix_A=zeros(N,N);

matrix_A(1,1)= 1 + lambda + lambda * h * A;//Robin BC
matrix_A(1,2)= -lambda;//Robin BC

matrix_A(N,N)= 1 + lambda;//Neumann BC
matrix_A(N,N-1)= -lambda;//Neumann BC

for i=2:N-1
    matrix_A(i,i-1)= -lambda/2;
    matrix_A(i,i)= 1 + lambda;
    matrix_A(i,i+1)= -lambda/2;
end

//MATRIX B
matrix_B=zeros(N,N);

matrix_B(1,1)=1 - lambda - lambda * h * A;
matrix_B(1,2)= lambda;

matrix_B(N,N)= 1 - lambda;
matrix_B(N,N-1)= lambda;

for i=2:N-1
    matrix_B(i,i-1)= lambda/2;
    matrix_B(i,i)= 1 - lambda;
    matrix_B(i,i+1)= lambda/2;
end


//SOLVE AND PLOT
C=(matrix_A^(-1))*matrix_B;
for i=1:length(TT)
    U=C*Up';
    Up=U';
    drawlater();    //not bad to get some continuous render...
    clf();
    a=gca();

    MM(i)=max(U);
    str = [ "time = "+string(arr(TT(i),3))+"      max = "+string(arr(MM(i),2))+"      max_pos = "+string(arr(find_max_pos(U),3)) ] ;
    xstring( max(XX)/5, u0max/2, str ) ;
    e = gce() ;
    e.font_angle = 0 ;
    e.font_size    = 4  ;
    e.font_style   = 0 ;
    //e.box = 'on' ;
    a.title.text="Heat eq Robin half-line (alpha = "+string(alpha)+") (Crank-Nicolson scheme)";
    a.title.font_size = 5;
    a.data_bounds = [0, -u0max/10 ; L, u0max]
    plot(XX',Up');
    xfarc(find_max_pos(U),0,max(XX)/80,max(u0max)/30,64*360,64*360)
    e=gce()
    e.background=4  //cyan
    f=gcf();
    f.figure_position = [-8,-8]
    f.figure_size = [1936,1096]
    drawnow();
    //xs2png(f,"save/t_"+string(i)+".png");

    Up=U';
end