function hdg
  omega=10;
  number_media=3;
  length_media=[1,1,1];
  mu_media=[1, 1, 1];
  rho_media=[1, 1, 1];
  point_per_wavelength=10;
  alpha=0.;
  point_per_disc=20;
%si source ponctuelle, position de la source
  x_source=12;
  
  [nbele,nbpoints,ele,coor,rho_cell,mu_cell]=create_mesh(omega,number_media,length_media,rho_media,mu_media,point_per_wavelength);
  order=1;
  for p=order
    basis(p)=basis_functions(p);
  end
  for i=1:nbele
    node=ele(i,:);
    x0=coor(node(1));
    x1=coor(node(2));
    rho=rho_cell(i);
    mu=mu_cell(i);
    [matrix(i).coeff, matrix(i).rhs]=create_matrix(x0,x1,rho,mu,order,alpha,basis(order),omega,x_source);
  end
  A=sparse(nbpoints,nbpoints);

 %Condition de bord gauche 
  A(1,1)=matrix(1).coeff(1,1)-2*alpha*matrix(1).coeff(p+2,1)+2*alpha;
   A(2,1)=matrix(1).coeff(1,2)-2*alpha*matrix(1).coeff(p+2,2);

    B(1)=-matrix(1).rhs(p+1)+2*alpha*(matrix(1).rhs(2*p+2));

   for i=2:nbpoints-1
    A(i-1,i)=-matrix(i-1).coeff(p+1,1)-alpha*matrix(i-1).coeff(2*p+2,1);
    A(i,i)=-matrix(i-1).coeff(p+1,2)-alpha*matrix(i-1).coeff(2*p+2,2)+matrix(i).coeff(1,1)-alpha*matrix(i).coeff(p+2,1)+2*alpha;
    A(i+1,i)=matrix(i).coeff(1,2)-alpha*matrix(i).coeff(p+2,2);
    B(i)=-(matrix(i).rhs(p+1)-matrix(i-1).rhs(1))+alpha*(matrix(i).rhs(2*p+2)+matrix(i-1).rhs(p+2));
  end
    B(nbpoints)=matrix(nbpoints-1).rhs(1)+2*alpha*matrix(nbpoints-1).rhs(p+2);

  % Condition de bord droit, 
   A(nbpoints-1,nbpoints)=-matrix(nbpoints-1).coeff(p+1,1)-2*alpha*matrix(nbpoints-1).coeff(2*p+2,1);
  A(nbpoints,nbpoints)=-matrix(nbpoints-1).coeff(p+1,2)-2*alpha*matrix(nbpoints-1).coeff(2*p+2,2)+2*alpha;
  % Condition de bord gauche et droit,  à commenter si source ponctuelle
  B(1)=0;
  B(nbpoints)=sin(omega*coor(nbpoints));
  lambda=B/A;
			 %  plot(imag(lambda));
 % figure
 %plot(real(lambda));
			 %  for i=1:nbpoints
 % x=coor(i);
  % lambda(i)=sin(omega*x);
 % end
  %plot(imag(A*lambda'))
 
  nb=0;
  for i=1:nbele
    node=ele(i,:);
    lambda_K(1,1)=lambda(node(1));
    lambda_K(2,1)=lambda(node(2));
    W_K= matrix(i).coeff*lambda_K+matrix(i).rhs;
    P_K=W_K(1:p+1);
    U_K=W_K(p+2:2*p+2);
    for j=1:point_per_disc+1
      nb=nb+1;
      x=(j-1)/point_per_disc;
      X(nb)=coor(node(1))+x*(coor(node(2))-coor(node(1)));
      P(nb)=0;
      U(nb)=0;
      for k=1:p+1
	      val=polyval(basis(p).phi(k).coeff,x);
	      P(nb)=P(nb)+P_K(k)*val;
	      U(nb)=U(nb)+U_K(k)*val;
      end
    end
  end
  plot(X,imag(P))
  hold on
  plot(X,imag(U))
