// G. Chiavassa juin 2011
//methode de tir lineaire
//equation T"= alpha T
// CL1: T(0) et T(1) fixees
// CL2: T(0) et T'(1) fixees


function DU=systeme(t,U),

DU = A*U;
endfunction;

function y=solexacte_CL1(x,alpha,T0,T1)

 mu=(T1-T0*exp(sqrt(alpha)))/(exp(-sqrt(alpha))-exp(sqrt(alpha)));
 lambda=T0-mu;

 y= lambda*exp(x*sqrt(alpha)) + mu* exp(-x*sqrt(alpha));

endfunction;

function y=solexacte_CL2(x,alpha,T0)

 mu=T0*exp(sqrt(alpha))/(exp(-sqrt(alpha))+exp(sqrt(alpha)));
 lambda=T0-mu;

 y= lambda*exp(x*sqrt(alpha)) + mu* exp(-x*sqrt(alpha));

endfunction;




//prog principal

N=200;
h=1./N;
x=[0:h:1];

alpha=10.;

A=[0 1; alpha 0];

/////CL1:  T(0)=TO et T(1)=T1; /////

T0=300;
T1=200;

// sol avec U0=[T0 0]';

 
 U0=[T0 0]'
 Y1=ode("rk",U0,0.,x,systeme);

 xset("window",1);
   clf();
   title(" solution CL1: Y1");
   plot2d(x,Y1(1,:));

// sol avec U0=[T0 1]';

 U0=[T0 1]'
 Y2=ode("rk",U0,0.,x,systeme);
 xset("window",2);
   clf();
   title(" solution CL1: Y2");
   plot2d(x,Y2(1,:));

// calcul des coeffs:

   c1=(T1-Y2(1,$))/(Y1(1,$)-Y2(1,$));
   c2=(Y1(1,$)-T1)/(Y1(1,$)-Y2(1,$));

   U=c1*Y1(1,:) +c2*Y2(1,:);
   U($)
   xset("window",3);
   clf();
   title(" solution CL1: U");
   plot2d(x,U,3);
   y=solexacte_CL1(x,alpha,T0,T1);
   plot2d(x,y,-1);


/////CL2:  T(0)=TO et T'(1)=0; /////


  b1= -Y2(2,$)/(Y1(2,$)-Y2(2,$));
  b2= Y1(2,$)/(Y1(2,$)-Y2(2,$));

  U=b1*Y1(1,:) +b2*Y2(1,:);
  
   xset("window",4);
   clf();
   title(" solution CL2: U");
   plot2d(x,U,3);
   z=solexacte_CL2(x,alpha,T0);
   plot2d(x,z,-1);





