function [ PHISOL ] = calculgrambis(Nbtermes,NPT,T,dt,V,z0,dz0,alpha0,dalpha0,beta0)
m=1;
ro=1.2;
J0=3;
k=1;
c=2;
a0=1.;
a=-.2;
S=1.;
L=.7;
bL=2.;
%========
% Construction de M K et F
%========
M=[m, a*cos(alpha0);a*cos(alpha0),J0];
K=[k ,-ro*S*V^2*dcz0(alpha0)/2;0 c-ro*S*L*dcm0(alpha0)/2];
Rad=linsolve(M,K');
F0=[ro*S*V^2*cz0(alpha0)/2,ro*S*L*V^2*cm0(alpha0)/2];
B0=[ro*S*V^2*dczg(beta0)/2,ro*S*L*V^2*dcmg(beta0)/2];
C0=[ro*S*V*(cxg(beta0)+dczg(beta0))/2,-bL*ro*S*L*V*(2*cmg(beta0)+dcmg(beta0))/2];

C=linsolve(M,C0');
B=linsolve(M,B0');
F=linsolve(M,F0');
%==================
% donnees initiales
%==================
% construction de la matrice de Gram
GRAMM=zeros(4,4);
%=======
phi0=linspace(1,4,4);
deltaphi0=linspace(1,4,4);
Lsec=linspace(1,2,4)*0.;
     P1dot=linspace(1,2,NPT);
     P2dot=linspace(1,2,NPT);
     Q1dot=linspace(1,2,NPT);
     Q2dot=linspace(1,2,NPT);
%===========
for k1=1:4
  for i1=1:4
        if(i1 == k1)
                phi0(i1)=1;
        else
                phi0(i1)=0;
        end
  end
    [P1,P2]=soluedo(Rad,NPT,dt,phi0);
     P1dot(1)=phi0(2);
     P2dot(1)=phi0(4);
     phii0=M*phi0;
     for i=2:NPT
         P1dot(i)=(P1(i)-P1(i-1))/dt;
         P2dot(i)=(P2(i)-P2(i-1))/dt;
     end
     Lsec(k1)=dz0*phii0(1)+dalpha0*phii0(3)-z0*phii0(2)-alpha0*phii0(4);
     squatro=0;
     for jkl=1:NPT
         squatro=squatro+F(1)*P1(jkl)+F(2)*P2(jkl);
     end
     Lsec(k1)=Lsec(k1)+squatro*dt;
  %==============
  for k2=k1:4
       for i2=1:4
          if(i2 == k2)
               deltaphi0(i2)=1;
          else
               deltaphi0(i2)=0;
          end
        end
%======
% Remplissage des solution elementaire
%======
    [Q1,Q2]=soluedo(Rad,NPT,dt,deltaphi0);
     Q1dot(1)=deltaphi0(2);
     Q2dot(1)=deltaphi0(4);
     for i=2:NPT   
      Q1dot(i)=(Q1(i)-Q1(i-1))/dt;
      Q2dot(i)=(Q2(i)-Q2(i-1))/dt;
     end
%========   
% calcul des termes de GRAM
%========
  
    
        coef=dt/(a0);
        sbis=0;
      
        for jj=1:NPT
      
        sbis=sbis+(C0(1)*P1dot(jj)+C0(2)*P2dot(jj)-B0(1)*P1(jj)-B0(2)*P2(jj))*(C0(1)*Q1dot(jj)+C0(2)*Q2dot(jj)-B(1)*Q1(jj)-B(2)*Q2(jj));
        end
        
   
    GRAMM(k1,k2)=GRAMM(k1,k2)+sbis*coef;
    GRAMM(k2,k1)=GRAMM(k1,k2);
 end
end
%======
%  calcul de PHI
%======
PHISOL=linsolve(GRAMM,Lsec');


end