messagebox(['Integrazione problema n.1  con' 'schemi alle differenze finite.' ' '  'Ver 3.1  (24/3/2010)'], 'Program3',"modal", "Ok" );


cont0="Y";
while cont0=="Y",

pi=atan (1)*4;

txt=['mu    '];
sig=x_mdialog('Parametri equazione  di diffusione',txt,['1.0']);
mu=evstr(sig(1));
lambda=0;
//mu=1;

N=51;
dt=0.0000499*4;
itermax=100;
Resmin=1e-5;

cambia=2;
while cambia==2, 
sN=msprintf('%i', N);
sdt=msprintf('%f', dt);
sitermax=msprintf('%i', itermax);
sResmin=msprintf('%e', Resmin);
txt=['N. nodi';'dt    ';'N. max iter';'Res. min'];
//sig=x_mdialog('Parametri  integrazione ',txt,[ 11 ;'0.004';'50';'1e-5']);
sig=x_mdialog('Parametri  integrazione ',txt,[ sN ; sdt ; sitermax ; sResmin]);
N=evstr(sig(1));
dt=evstr(sig(2));
itermax=evstr(sig(3));
Resmin=evstr(sig(4));

//N=101;
dx=1/(N-1);
c=lambda*dt/dx;
s=mu*dt/dx^2;


tschema=['FTCS (esplicito,  I ord. tempo, II ord. spazio)';...
         'BTCS (implicito,  I ord. tempo, II ord. spazio)';...
         'CN   (implicito, II ord. tempo, II ord. spazio)'];
nschema=x_choose_modeless(tschema,['Scegli lo schema alle differenze finite'])


testo(1) =strcat( ['Riassunto dei parametri equazione ed  integrazione:']);
testo(2)=strcat( ['------------------------------------------------------------']);
testo(3) =strcat( ['Parametri equazione ----->  lambda:' string(msprintf('%2.5f',lambda   ))']);
testo(4) =strcat( ['                                mu:' string(msprintf('%2.5f',mu       ))']);
testo(5) =strcat( ['Parametri integrazione --> N. nodi:' string(msprintf('%2.5f',N   ))']);
testo(6) =strcat( ['                                dx:' string(msprintf('%2.5f',dx  ))']);
testo(7) =strcat( ['                                dt:' string(msprintf('%2.5f',dt  ))']);
testo(8) =strcat( ['                                 c:' string(msprintf('%2.5f',c   ))']);
testo(9) =strcat( ['                                 s:' string(msprintf('%2.5f',s   ))']);
testo(10) =strcat( ['                       N. iter max:' string(msprintf('%2.5f',itermax))']);
testo(11)=strcat( ['                         Res. min.:' string(msprintf('%2.5f',Resmin))']);
testo(12)=strcat( ['            Schema di integrazione:' string(msprintf('%s',tschema(nschema) ))']);
testo(13)=strcat( ['------------------------------------------------------------']);
testo(14)=strcat( ['Vuoi cambiare i parametri?']);
cambia= messagebox(testo,'Program3',"question",['No' 'Yes'], "modal");
end,


X=[0 : dx : 1.]';
gn=sin(pi*X);
gnp1=zeros(1,N)';

iter=2000;


if nschema==1 then
  //FTCS
  A=zeros(N,N);
  A(1,1)=1.;
  for i=2:N-1
    A(i,i+1)=s-c/2;
    A(i,i)=1-2*s;
    A(i,i-1)=(c/2+s);
  end
  A(N,N)=1.;
  L=A;
elseif nschema==2 then
  //BTCS
  A=zeros(N,N);
  A(1,1)=1.;
  for i=2:N-1
    A(i,i+1)=c/2-s;
    A(i,i)=1+2*s;
    A(i,i-1)=-(c/2+s);
  end
  A(N,N)=1.;
  L=inv(A);
else
  //CN
  A=zeros(N,N);
  B=zeros(N,N);
  A(1,1)=1.;
  B(1,1)=1.;
  for i=2:N-1
    A(i,i+1)=1/2*(c/2-s);
    B(i,i+1)=1/2*(s-c/2);
    A(i,i)=1+s;
    B(i,i)=1-s;
    A(i,i-1)=-1/2*(c/2+s);
    B(i,i-1)=1/2*(c/2+s);
  end
  A(N,N)=1.;
  B(N,N)=1.;
  L=inv(A)*B;
end,



win0=scf(0);
clf(win0,"reset");
xtitle('Evoluzione della soluzione', 'x', 'g');
xset("window",0)
xselect();

xselect();
plot2d(X,gn,-1,rect=[0.,0.,1.,1.1]);
s1=gce();



//xset("mark",1,2);
t=0.;

ge=sin(pi*X)*exp(-pi^2*mu*t);
plot2d(X,ge,2,rect=[0.,0.,1.,1.1]);
s2=gce();

res=1e+39;
gnp1max=0.;
nn=1;
while nn <= itermax & res > Resmin & gnp1max <= 1.05 , 
  gnp1=L*gn;
  gnp1max=max(gnp1);
  if gnp1max > 1.5 then
   messagebox(['Attenzione   g>1','Probabile schema di integrazione instabile'], 'Program3',"modal", "Ok" );
   break;
  end;

    s1.children.data=[X gnp1];
  
  t=t+dt;
  ge=sin(pi*X)*exp(-pi^2*mu*t);
  s2.children.data=[X ge];

  gn=gnp1;
  nn=nn+1;
end
xselect();



cont1=x_dialog('Lettura soluzione in un punto?','Y');
if cont1=="Y" then
 cont2="Y";
 while cont2=="Y",
 inode=evstr(x_dialog('Scegli il nodo da valutare(1-N)','0'));
 testo1(1) =strcat( ['Confronto tra la soluzione numerica e la soluzione esatta in']);
 testo1(2) =strcat( ['                    t: ' string(msprintf('%f', t           ))']);
 testo1(3) =strcat( ['                    X: ' string(msprintf('%f', X(inode)    ))']);
 testo1(4) =strcat( ['----------------------------------------------------------------']);      
 testo1(5) =strcat( ['soluzione -->  Esatta:' string(msprintf('%f',ge(inode)     ))']);
 testo1(6) =strcat( ['             Numerica:' string(msprintf('%f',gnp1(inode)   ))']);
 testo1(7) =strcat( ['----------------------------------------------------------------']);      
 testo1(8) =strcat( ['Soluzione numerica ottenuta con:']);
 testo1(9) =strcat( [ string(msprintf('%s',tschema(nschema))) ]);
 testo1(10) =strcat( ['integrando con  dx:'  string(msprintf('%f',dx )) '      dt:' string(msprintf('%f',dt ))]);
 messagebox(testo1,'Risultato',"modal", "Ok");
 cont2=x_dialog('Valutazione di un altro nodo ?','Y');
 end,
end,
cont0=x_dialog('Nuova valutazione ?','Y');
end,



quit;