%   asiatiques.m                                 Francois Dubois  
%   Options asiatiques  avec Octave.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
c0 =  sprintf (' "asiatiques"  exits also in language C ; ask to Francois Dubois !' )
%     asia1     test de la derivee seconde              28, 30  mars 2003
%     asia2     premieres experiences              30 mars, 1 avril  2003
%     asia3     test des nouvelles derivees secondes     7, 8 avril  2003
%     asia4     schema d ordre 3 en espace                  8 avril  2003
%     asia5     post-traitement soigne en x=1   9, 14, 18, 21 avril  2003
%     asia6     schema d ordre 4 en temps          21, 25, 28 avril  2003
%     asia7     changement regularisant     28, 29 avril, 01, 03 mai 2003
%     asia8     presentation unique pour l ordre 2           04  mai 2003
%     asia9     version web                                  09  mai 2003
%     asiatiques    correctins mineures                     22 avril 2006
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   ordre en temps      2    schema de Crank Nicolson
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   choix de l approche
%   app = 1         approche initiale  par differences finies 
%   app = 2         changement de variables et masse de Dirac
    app = 2
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%   taux d interet  r  et volatilite     sigma
%                  0.02                   0.30          Shi  facile 
%                  0.02                   0.05          Shi  difficile 
%                  0.15                   0.05          Zvan et al   
%                  0.10                   0.20          Lapeyre-Temam   
    r     = 0.02
    sigma = 0.05
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
%   nombre de pas de temps                      nt
%   nombre de sommets a l instant initial       nj
%   nombre total de sommets                     ns=nt + nj
%   temps final                                 tf   
    nt    = 100
    nj    = nt
    ns    = nt + nj  ;           
    tf    = 1      ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%                                                          ecritures
    [chordre,msg]  = sprintf  (' approach  : %1i', app) ; 
    [chr,msg]   = sprintf  ('    rate  : %6.3f',r) ; 
    [chsig,msg] = sprintf  ('    volatility  : %5.2f',sigma) ; 
    [chnt,msg]  = sprintf  ('    nt  : %4i',nt) ;
    [chntg,msg]  = sprintf  ('   nb of time steps  : %4i',nt) ;
    [chnj,msg]  = sprintf  ('    nj  : %4i',nj) ;
    [chtf,msg]  = sprintf  ('    final time  : %5.2f',tf);
    titx = [chordre,chr , chsig , chntg , chnj] ; titx ; 
    titx1 = [chordre,chr , chsig , chtf] ; titx1 ; 
    title(c0) ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%   calcul du pas de temps dt     
    dt    = tf / nt
    'tap return to see the results around the value "1"'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%   grille geometrique en y 
%   ro   = exponentielle de r multiplie par dt
%   y0   = exponentielle de -r multiplie par tf
    ro   = exp(r*dt)  ; 
    rom1 = ro -1    ;
    y0   = exp(-r*tf)  ;
    y  = zeros(ns,1) ;
    y(1) = y0 * ro ;
    for j = 2:ns ;  y(j) = ro * y(j-1) ;  end ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                           derivee seconde centree
    y21  = 1 / (ro*ro) ; 
    y31  = (1+ro) / (ro*ro) ; 
    y41  = (1+ro+ro*ro) / (ro*ro) ; 
    y51  = (1+ro)*(1+ro*ro) / (ro*ro) ;
    y32  = 1 / ro ;
    y42  = (1+ro) / ro ;  
    y52  = (1+ro+ro*ro) / ro ;
    y43  = 1 ;   
    y53  = 1+ro ;   
    y54  = ro ;  
    yst1  = -y31 ;
    yst2  = -y32 ; 
    yst4  =  y43 ; 
    yst5  =  y53 ; 
%
    d2c    = zeros(5,1) ;
    d2c(1) =  2 * (yst4*yst5 + yst2*(yst4+yst5)) / (y21*y31*y41*y51) ; 
    d2c(2) = -2 * (yst4*yst5 + yst1*(yst4+yst5)) / (y21*y32*y42*y52) ;
    d2c(3) =  2 * (yst1*yst2 + yst1*yst4 + yst1*yst5 + yst2*yst4 \
                  + yst2*yst5 + yst4*yst5 ) / (y31*y32*y43*y53) ;
    d2c(4) = -2 * (yst2*yst5 + yst1*(yst2+yst5)) / (y41*y42*y43*y54) ;
    d2c(5) =  2 * (yst2*yst4 + yst1*(yst2+yst4)) / (y51*y52*y53*y54) ; 
    for  i = 1:5 ; d2c(i) = d2c(i) / ((ro-1)*(ro-1)) ; end ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%                              derivee seconde decentree vers la gauche  
    y21  = 1 / ro ; 
    y31  = (1+ro) / ro ; 
    y41  = (1+ro+ro*ro) / ro  ; 
    y51  = (1+ro)*(1+ro*ro) / ro ;
    y32  = 1  ;
    y42  = (1+ro)  ;  
    y52  = (1+ro+ro*ro)  ;
    y43  = ro ;   
    y53  = ro + ro*ro ;   
    y54  = ro * ro  ;  
    yst1  = -y21 ;
    yst3  =  y32 ; 
    yst4  =  y42 ; 
    yst5  =  y52 ; 
%
    d2dg    = zeros(5,1) ;
    d2dg(1) =  2 * (yst3*yst4 + yst4*yst5 + yst5*yst3) / (y21*y31*y41*y51) ; 
    d2dg(2) = -2 * (yst1*yst3 + yst1*yst4 + yst1*yst5 + yst3*yst4 \
                  + yst3*yst5 + yst4*yst5 ) / (y21*y32*y42*y52) ;
    d2dg(3) =  2 * (yst1*yst4 + yst4*yst5 + yst5*yst1) / (y31*y32*y43*y53) ;
    d2dg(4) = -2 * (yst1*yst3 + yst3*yst5 + yst5*yst1) / (y41*y42*y43*y54) ;
    d2dg(5) =  2 * (yst3*yst4 + yst4*yst1 + yst1*yst3) / (y51*y52*y53*y54) ;
    for  i = 1:5 ; d2dg(i) = d2dg(i) / ((ro-1)*(ro-1)) ; end ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%                                      derivee seconde au bout a droite  
%          attention    avec une condition de Neumann homogene a droite
    y21  =  1 / (ro*ro*ro) ; 
    y31  = (1+ro) / (ro*ro*ro) ; 
    y41  = (1+ro+ro*ro) / (ro*ro*ro)  ; 
    y32  =  1 / (ro*ro) ;
    y42  = (1+ro) / (ro*ro)  ;  
    y43  = 1 / ro ;   
    yst1  = -y41 ;
    yst2  = -y42 ; 
    yst3  = -y43 ; 
%
    d2pd    = zeros(4,1) ;
    d2pd(1) =  2 * (y42*y43) / (y21*y31*y41*y41) ; 
    d2pd(2) = -2 * (y41*y43) / (y21*y32*y42*y42) ;
    d2pd(3) =  2 * (y41*y42) / (y31*y32*y43*y43) ;
    d2pd(4) = -2 * (yst1*yst1*yst2*yst2 + yst2*yst2*yst3*yst3         \ 
            + yst3*yst3*yst1*yst1 + yst1*yst2*yst3*(yst1+yst2+yst3)) \ 
                                    / (y41*y41*y42*y42*y43*y43) ;
     for  i = 1:4 ; d2pd(i) = d2pd(i) / ((ro-1)*(ro-1)) ; end ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%              derivee seconde a l avant dernier point au bout a droite  
%           attention:  avec une condition de Neumann homogene a droite
    y21  =  1 / (ro*ro) ; 
    y31  = (1+ro) / (ro*ro) ; 
    y41  = (1+ro+ro*ro) / (ro*ro)  ; 
    y32  =  1 / (ro) ;
    y42  = (1+ro) / (ro)  ;  
    y43  = 1  ;   
    yst1  = -y41 ;
    yst2  = -y42 ; 
    yst3  = -y43 ; 
%
    d2pad    = zeros(4,1) ;
    d2pad(1) =  2 * (3*yst3 - 2*yst2)*yst3 / (yst1*yst1*y21*y31) ; 
    d2pad(2) = -2 * (3*yst3 - 2*yst1)*yst3 / (yst2*yst2*y21*y32) ; 
    d2pad(3) =  2 * (6*yst3*yst3 - 3*yst3*(yst1+yst2) + yst1*yst2 ) \
                            /  (yst3*yst3*y32*y31) ; 
    d2pad(4) = -2 * (3*yst3*yst3*yst3*(yst1+yst2) +                 \
                     yst1*yst2-2*(yst1*yst1+yst2*yst2)*yst3*yst3   \
               -2*yst1*yst2*(yst1+yst2)*yst3 + yst1*yst1*yst2*yst2) \
                           /  (yst1*yst1*yst2*yst2*yst3*yst3) ; 
    for  i = 1:4 ; d2pad(i) = d2pad(i) / ((ro-1)*(ro-1)) ; end ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                    condition initiale
    uu   = zeros(ns,1) ;
    f0   = zeros(ns,1) ;
    sm   = zeros(ns,1) ;
    zz   = zeros(ns,1) ;
    zp   = zeros(ns,1) ;
    xxx  = zeros(ns,1) ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                     condition limite et masse de Dirac
    cl       = zeros(nt+1,1) ; 
    delta    = zeros(nt,1) ; 
%    
    if (app==1) ;
    for    k = 1:nt ; 
    cl(k)    = (1-y(nt-k+1)) / (r*tf) ;   delta(k) = 0 ; end ; 
    cl(nt+1) =  (1-y0) / (r*tf) ; 
    else ; 
    for    k = 1:nt ; 
    cl(k)    = 0  ; 
    time     = (k-1)*dt ; 
    delta(k) = sigma*sigma*(dt +2*exp(-r*time)*(exp(-r*dt)-1)/r  \
                       - exp(-2*r*time)*(exp(-2*r*dt)-1)/ (2*r)) \
               / (r*tf*(y(nt+1)-y(nt-1)))  ;   end ;   
    end ;  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                           iterations en temps 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  
    tt      = zeros(nt+1,1) ;
    time    = 0 ;
    extps   = exp(-r*time) ; 
    for   j = 1 : ns
    zz(j)   = sigma*sigma* (y(j)-extps)*(y(j)-extps) / (2*y(j)*y(j)) ; 
    zp(j)   = r*sigma*sigma*(y(j)-extps) /(y(j)*y(j)) ;     end ;
%
%                                          k     numero du pas de temps 
%                j0 :  adresse du point numero zero a prendre en compte
%                           nddl : ordre du systeme lineaire a resoudre
    for   k = 1:nt ;
    tt(k)   = time ; 
    extps   = exp(-r*time) ; 
    extpsp1 = exp(-r*time) * exp(-r*dt) ; 
    j0      = nt - k  ;
    nddl    = ns - j0  ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                          calcul de f0 
%
%                                         au premier point, tout est nul
    j  = j0+1 ; f0(j) = 0 ; 
%                                                           second point
    j     = j0+2 ; 
    f0(j) = ( d2dg(1)*cl(k) + d2dg(2)*uu(j) + d2dg(3)*uu(j+1) + \
              d2dg(4)*uu(j+2) ) * zz(j) ;
%
%                                                          point courant
    for l = 3:nddl-2 ;  j = j0+l ; 
    f0(j) = ( d2c(1)*uu(j-2) +  d2c(2)*uu(j-1) + d2c(3)*uu(j) +   \
              d2c(4)*uu(j+1) +  d2c(5)*uu(j+2) ) * zz(j) ; end ; 
%
%                           derivee seconde speciale pour le point nddl-1
    l     = nddl-1 ; j = j0+l ; 
    f0(j) = ( d2pad(1)*uu(j-2) + d2pad(2)*uu(j-1) + d2pad(3)*uu(j) +   \
              d2pad(4)*uu(j+1) ) * zz(j) ;
%
%                          derivee seconde speciale pour le dernier point 
    l = nddl ; j = j0+l ; 
    f0(j) = ( d2pd(1)*uu(j-3) + d2pd(2)*uu(j-2) + d2pd(3)*uu(j-1) +  \
              d2pd(4)*uu(j) ) * zz(j) ;
%  f0
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                       second membre pour le schema de Crank Nicolson 
    for   j = 1 : ns
    zz(j)   = sigma*sigma*(y(j)-extpsp1)*(y(j)-extpsp1) / (2*y(j)*y(j)) ; 
    zp(j)   = r*sigma*sigma*(y(j)-extpsp1) /(y(j)*y(j)) ;     end ;
    j0      = nt - k  ;
    for l = 1:nddl ; j = j0+l ;  
    sm(l) = uu(j) + dt*f0(j)/2 ; end ;  
    j1    = j0+1 ;  sm(1) = sm(1) + dt * d2dg(1)*cl(k+1)*zz(j1) / 2 ;
    j2    = j0+2 ;  sm(1) = sm(1) + dt * d2c(1) *cl(k+1)*zz(j2) / 2 ;
%
%                                 ajout de la masse de Dirac en y = 1
    j = nt ; l = j - j0 ; 
    sm(l) = sm(l) + delta(k) ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                          matrice quasi-pentadiagonale 
    aa = zeros(nddl,1) ;  bb = zeros(nddl,1) ; cc = zeros(nddl,1) ;      
    dd = zeros(nddl,1) ;  ee = zeros(nddl,1) ; 
%                                                         premier point
    j1   = j0+1       ; 
    aa(1) = 1 -  d2dg(2)*dt*zz(j1) / 2  ; 
    bb(2) =   -  d2dg(3)*dt*zz(j1) / 2  ; 
    dd(3) =   -  d2dg(4)*dt*zz(j1) / 2  ;  
    a14   =   -  d2dg(5)*dt*zz(j1) / 2  ; 
%                                                         second  point
    j2   = j0+2       ; 
    cc(2) =   -  d2c(2)*dt*zz(j2) / 2  ;
    aa(2) = 1 -  d2c(3)*dt*zz(j2) / 2  ;
    bb(3) =   -  d2c(4)*dt*zz(j2) / 2  ;
    dd(4) =   -  d2c(5)*dt*zz(j2) / 2  ;
%                                                   points 3 a nddl - 2   
    for l = 3:nddl-2 ;   j = j0 + l ;     
    ee(l)   =   -    d2c(1)*dt*zz(j) / 2  ;
    cc(l)   =   -    d2c(2)*dt*zz(j) / 2  ;
    aa(l)   = 1 -    d2c(3)*dt*zz(j) / 2  ;
    bb(l+1) =   -    d2c(4)*dt*zz(j) / 2  ;
    dd(l+2) =   -    d2c(5)*dt*zz(j) / 2  ;
    end ;
%                                                        point  nddl - 1   
    l = nddl-1 ;   j = j0 + l ;      
    ee(l)   =   -    d2pad(1)*dt*zz(j) / 2  ;
    cc(l)   =   -    d2pad(2)*dt*zz(j) / 2  ;
    aa(l)   = 1 -    d2pad(3)*dt*zz(j) / 2  ;
    bb(l+1) =   -    d2pad(4)*dt*zz(j) / 2  ;
%
%                                                            point  nddl    
    l = nddl ;   j = j0 + l ;      
    a41     =   -    d2pd(1)*dt*zz(j) / 2  ;
    ee(l)   =   -    d2pd(2)*dt*zz(j) / 2  ;
    cc(l)   =   -    d2pd(3)*dt*zz(j) / 2  ;
    aa(l)   = 1 -    d2pd(4)*dt*zz(j) / 2  ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                        factorisation l*u de la  matrice pentadiagonale  
    aal   = zeros(nddl,1) ;            
    bbu   = zeros(nddl,1) ;           
    ccl   = zeros(nddl,1) ;          
    ddu   = zeros(nddl,1) ;          
    eel   = zeros(nddl,1) ;  
%
    aal(1)  = aa(1) ; 
    bbu(2)  = bb(2) / aal(1) ;
    ccl(2)  = cc(2) ; 
    aal(2)  = aa(2) - ccl(2) * bbu(2) ;   
    uu14    = a14 / aal(1) ; 
    for  i = 3:nddl ;
    eel(i)  = ee(i) ;
    ddu(i)  = dd(i) / aal(i-2) ;
    ccl(i)  = cc(i) - eel(i)*bbu(i-1) ; 
    bbu(i)  = (bb(i) - ccl(i-1)*ddu(i)) / aal(i-1) ;
    aal(i)  = aa(i) - ddu(i)*eel(i) - bbu(i)*ccl(i) ;
    end ;
    ll41    = a41 ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                    descente-remontee        
    xx  = zeros(nddl,1) ;  yy = zeros(nddl,1) ; 
%          
    yy(1)  = sm(1) / aal(1) ;
    yy(2)  = (sm(2) - ccl(2)*yy(1)) / aal(2) ;  
    for  i = 3:nddl-1 ;
    yy(i)  = (sm(i) - ccl(i)*yy(i-1) - eel(i)*yy(i-2)) / aal(i) ; end ;
    i      = nddl ; 
    yy(i)  = (sm(i) - ccl(i)*yy(i-1) - eel(i)*yy(i-2) - ll41*yy(i-3)) \
                              / aal(i) ;
%
    xx(nddl)   = yy(nddl) ;
    xx(nddl-1) = yy(nddl-1) - bbu(nddl)*xx(nddl) ;  
    for      i = 3:nddl ; j = nddl+1-i ;
    xx(j)      = yy(j) - bbu(j+1)*xx(j+1) - ddu(j+2)*xx(j+2) ; end ; 
    xx(1)      = xx(1) - uu14*xx(4) ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%                               nouvelle valeur au nouveau pas de temps 
%
    vv     = zeros(ns,1) ;
    sol    = zeros(ns,1) ;
    phi    = zeros(ns,1) ;
%                                            ajout de ce qu on a enleve
    for  j = 1:nt ; 
    vv(j)  = (1-y(j)) / (r*tf) ; end 
%
    if (app==1) ;  
    for  j = 1:ns ;  uu(j)=cl(k+1) ; end ; end ; 
%
    for  i = 1:nddl ;  j = j0 + i ;  uu(j)  =   xx(i) ;   end ;
%
    if (app==1) ; 
    for  j=1:ns ;  phi(j) = uu(j) ; end ; 
    else ; 
    for  i = 1:nddl ;  j = j0 + i ;  phi(j) = xx(i) + vv(j) ; end ; end ; 
%
    for   i = 1:ns ; 
    xxx(i)  = ( y(i) * exp(r*(time+dt)) - 1) / (r*tf) ; end ; 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%     afffichage     %%%%%%%%%%%%%%%%%%%%%%%%
%
    [chdt,msg] = sprintf ('      dt numero : %3i',k);
    titx2 = [chr , chsig , chnt , chnj , chdt] ; titx2 ; 
    xlabel ( titx2 ) ;
    axis ([ 0.   2.05   -1  101 ]) ; 
%     gset xrange[0.: 2.05 ] ; gset yrange[-1 :   101 ] ;
%
    if (app==1) ; 
    plot (xxx,100*phi,'o3;schema initial;') 
    else ; 
    plot (xxx,100*phi,'o3;schema Dirac;') 
    end ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                            fin de la boucle en temps
    time = time + dt ; end ; 
    pause
    tt(k+1) = time ;   
    if (app==1) ; 
    for   i = 1:nddl ;  j = j0 + i ;  sol(j)  =  xx(i) ; end ; else ; 
    for   i = 1:nddl ;  j = j0 + i ;  
    sol(j)  =  vv(j) +  xx(i) ; end ; end ;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%     afffichage     %%%%%%%%%%%%%%%%%%%%%%%%
%    gset xrange[0.98: 1.05 ] ; gset yrange[-0.05:   3.05 ] ;
    axis ([ 0.98  1.05   -0.05    3.05  ]) ; 
    xlabel ( titx1 ) ;
     plot (xxx,100*sol,'x6;solution  calculee;'  )
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%     test = eps  % precision des calculs 
%                                               interpolation en x=1 
     yref = (1+r*tf)*y0  ; 
     zz   = log(1+r*tf)/log(ro) ; 
     j    = floor(zz) ;
     x0   = xxx(j)    ;
     x1   = xxx(j+1)  ;
     x10  = x1-x0     ;
     xi   = (1 - x0)/(x10);
     u0   = 100*sol(j)    ;
     u1   = 100*sol(j+1)  ;
     valun = (1-xi)*u0 + xi*u1 ; 
     [p1,msg] = sprintf ('     valeur interpolee affine =  %15.12f',valun) ;
     [p1,msg] = sprintf ('    affine interpolated value =  %15.12f',valun) ;
     p1
%                                                   interpolation p2
     x2   = xxx(j+2)  ;
     u2   = 100*sol(j+2)  ;
     x20  = x2-x0     ;
     x21  = x2-x1     ;
     a0   = (1-xi)*(x20-x10*xi)/(x20) ;
     b0   = xi*(x20-x10*xi)/(x21) ;
     c0   = -xi*(1-xi)*x10*x10/(x20*x21) ;
     val2g = a0*u0 + b0*u1 + c0*u2 ; 
     [p2,msg] = sprintf ('  valeur interpolee p2 gauche =  %15.12f',val2g) ;
     [p2,msg] = sprintf (' left   p2 interpolated value =  %15.12f',val2g) ;
     p2
%
     x2   = xxx(j-1)  ;
     u2   = 100*sol(j-1)  ;
     x20  = x2-x0     ;
     x21  = x2-x1     ;
     a0   = (1-xi)*(x20-x10*xi)/(x20) ;
     b0   = xi*(x20-x10*xi)/(x21) ;
     c0   = -xi*(1-xi)*x10*x10/(x20*x21) ;
     val2d = a0*u0 + b0*u1 + c0*u2  ;
     [p2,msg] = sprintf ('  valeur interpolee p2 droite =  %15.12f',val2d) ;
     [p2,msg] = sprintf (' right  p2 interpolated value =  %15.12f',val2d) ;
     p2
%                                                   interpolation p3
     x2  = xxx(j+2) ;     x3  = xxx(j+3)  ;
     u2  = 100*sol(j+2)  ; u3  = 100*sol(j+3)  ;
     x20 = x2-x0 ; x21 = x2-x1 ; x30 = x3-x0 ; x31 = x3-x1 ;x32 = x3-x2 ;
     a0   = (1-xi)*(x20-x10*xi)*(x30-x10*xi)/(x20*x30) ;
     b0   = xi*(x20-x10*xi)*(x30-x10*xi)/(x21*x31) ;
     c0   = -xi*(1-xi)*(x30-x10*xi)*x10*x10/(x20*x21*x32) ;
     d0   = xi*(1-xi)*(x20-x10*xi)*x10*x10/(x30*x31*x32) ;
     tt0 = a0+b0+c0+d0 ;
     tt1 = b0*x10 + c0*x20 + d0*x30 ; 1-x0 ;
     tt2 = b0*x10*x10 + c0*x20*x20 + d0*x30*x30 ; (1-x0)*(1-x0) ;
     tt3 = b0*x10*x10*x10 + c0*x20*x20*x20 + d0*x30*x30*x30 ;
     (1-x0)*(1-x0)*(1-x0)  ;
     val3gg = a0*u0 + b0*u1 + c0*u2 +d0*u3 ;
     [p3,msg] = sprintf ('  valeur interpolee p3 gauche =  %15.12f',val3gg) ;
     [p3,msg] = sprintf (' left   p3 interpolated value =  %15.12f',val3gg) ;
     p3
%                                                       p3 centree
     x2  = xxx(j-1) ;     x3  = xxx(j+2)  ;
     u2  = 100*sol(j-1)  ; u3  = 100*sol(j+2)  ;
     x20 = x2-x0 ; x21 = x2-x1 ; x30 = x3-x0 ; x31 = x3-x1 ;x32 = x3-x2 ;
     a0   = (1-xi)*(x20-x10*xi)*(x30-x10*xi)/(x20*x30) ;
     b0   = xi*(x20-x10*xi)*(x30-x10*xi)/(x21*x31) ;
     c0   = -xi*(1-xi)*(x30-x10*xi)*x10*x10/(x20*x21*x32) ;
     d0   = xi*(1-xi)*(x20-x10*xi)*x10*x10/(x30*x31*x32) ;
     val3cc = a0*u0 + b0*u1 + c0*u2 +d0*u3 ;
     [p3,msg] = sprintf ('  valeur interpolee p3 centre =  %15.12f',val3cc) ; 
     [p3,msg] = sprintf (' center p3 interpolated value =  %15.12f',val3cc) ;
     p3
%                                                       p3 droite
     x2  = xxx(j-2) ;     x3  = xxx(j-1)  ;
     u2  = 100*sol(j-2)  ; u3  = 100*sol(j-1)  ;
     x20 = x2-x0 ; x21 = x2-x1 ; x30 = x3-x0 ; x31 = x3-x1 ;x32 = x3-x2 ;
     a0   = (1-xi)*(x20-x10*xi)*(x30-x10*xi)/(x20*x30) ;
     b0   = xi*(x20-x10*xi)*(x30-x10*xi)/(x21*x31) ;
     c0   = -xi*(1-xi)*(x30-x10*xi)*x10*x10/(x20*x21*x32) ;
     d0   = xi*(1-xi)*(x20-x10*xi)*x10*x10/(x30*x31*x32) ;
     val3dd = a0*u0 + b0*u1 + c0*u2 +d0*u3 ;
     [p3,msg] = sprintf ('  valeur interpolee p3 droite =  %15.12f',val3dd) ;
     [p3,msg] = sprintf (' right  p3 interpolated value =  %15.12f',val3dd) ;
     p3
%                                                             p4 gauche 
     x0  = xxx(j)      ;  x1  = xxx(j+1)  ;
     x2  = xxx(j+2)    ;  x3  = xxx(j+3)   ;      x4  = xxx(j-1)  ;
     u2  = 100*sol(j+2)  ; u3  = 100*sol(j+3)  ; u4  = 100*sol(j-1)  ;
     x10 = x1-x0  ; x20 = x2-x0 ; x30 = x3-x0 ;  x40 = x4-x0 ;
     x21 = x2-x1  ; x31 = x3-x1 ; x41 = x4-x1 ;  x32 = x3-x2 ;  
     x42 = x4-x2  ; x43 = x4-x3  ; 
     a0  = (1-xi)*(x20-x10*xi)*(x30-x10*xi)*(x40-x10*xi) /(x20*x30*x40); 
     b0  = xi*(x20-x10*xi)*(x30-x10*xi)*(x40-x10*xi)  / (x21*x31*x41) ;
     c0  = -xi*(1-xi)*(x30-x10*xi)*(x40-x10*xi)*x10*x10 /(x20*x21*x32*x42); 
     d0  =  xi*(1-xi)*(x20-x10*xi)*(x40-x10*xi)*x10*x10 /(x30*x31*x32*x43); 
     e0  = -xi*(1-xi)*(x20-x10*xi)*(x30-x10*xi)*x10*x10 /(x40*x41*x42*x43); 
     val4gg = a0*u0 + b0*u1 + c0*u2 + d0*u3 + e0*u4 ;
     [p4,msg] = sprintf ('  valeur interpolee p4 gauche =  %15.12f',val4gg) ; 
     [p4,msg] = sprintf (' left   p4 interpolated value =  %15.12f',val4gg) ;
     p4
%                                                             p4 droite 
     x0  = xxx(j)      ;  x1  = xxx(j+1)  ;
     x2  = xxx(j+2)    ;  x3  = xxx(j-1)   ;      x4  = xxx(j-2)  ;
     u2  = 100*sol(j+2)  ; u3  = 100*sol(j-1)  ; u4  = 100*sol(j-2)  ;
     x10 = x1-x0  ; x20 = x2-x0 ; x30 = x3-x0 ;  x40 = x4-x0 ;
     x21 = x2-x1  ; x31 = x3-x1 ; x41 = x4-x1 ;  x32 = x3-x2 ;  
     x42 = x4-x2  ; x43 = x4-x3  ; 
     a0  = (1-xi)*(x20-x10*xi)*(x30-x10*xi)*(x40-x10*xi) /(x20*x30*x40); 
     b0  = xi*(x20-x10*xi)*(x30-x10*xi)*(x40-x10*xi)  / (x21*x31*x41) ;
     c0  = -xi*(1-xi)*(x30-x10*xi)*(x40-x10*xi)*x10*x10 /(x20*x21*x32*x42); 
     d0  =  xi*(1-xi)*(x20-x10*xi)*(x40-x10*xi)*x10*x10 /(x30*x31*x32*x43); 
     e0  = -xi*(1-xi)*(x20-x10*xi)*(x30-x10*xi)*x10*x10 /(x40*x41*x42*x43); 
     tt0 = a0+b0+c0+d0+e0 ;
     tt1 = b0*x10 + c0*x20 + d0*x30 + e0*x40 ;  1-x0 ;
     tt2 = b0*x10*x10 + c0*x20*x20 + d0*x30*x30  + e0*x40*x40 ;(1-x0)*(1-x0); 
     tt3 = b0*x10*x10*x10 + c0*x20*x20*x20 + d0*x30*x30*x30 + e0*x40*x40*x40; 
     (1-x0)*(1-x0)*(1-x0)  ;
     tt3 = b0*x10*x10*x10*x10 + c0*x20*x20*x20*x20 \ 
                              + d0*x30*x30*x30*x30 + e0*x40*x40*x40*x40 ;
     (1-x0)*(1-x0)*(1-x0)*(1-x0) ;
     val4dd = a0*u0 + b0*u1 + c0*u2 + d0*u3 + e0*u4 ;
     [p4,msg] = sprintf ('  valeur interpolee p4 droite =  %15.12f',val4dd) ; 
     [p4,msg] = sprintf (' right  p4 interpolated value =  %15.12f',val4dd) ;
     p4
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%






