function yy=varvec(y,r,x);
% -  March 11 2003(RMM)
% - Solves:    Dny+ ...+p Dy +q y = r(x)
% -  y=[y1 y2 ... yn]; are known solutions(linearly independent.
% -  A matrix form is used: A A1 A2 w=|A| W1=|A1| w2=|A2|  wn=|An|
% -  the method of variation is used to evaluate yp.
% - usage:
% - syms x r y ; 
% - y(1)=1;y(2)=x;y(3)=x^2;y(4)=x^3;r=1;
% - y=varvec(y,r,x)
% -
syms  w1 w u A0 A B yh  yp yg nsol yy 
nsol=size(y,2)
  for i=1:nsol
    disp(y(i));
  end;   
  disp(r)

   for j=1:nsol
      A0(1,j)=y(j);
      u(i,1)=0;
   end;
   for i=2:nsol
      for j=1:nsol
         A0(i,j)=diff(y(j),x,i-1);
      end;
   end;
   
   w=det(simplify(A0));
              
   for k=1:nsol
           for i=1:nsol
                 for j=1:nsol
                     B(i,j)=A0(i,j);
                    if j==k    B(i,j)=0;end;
          
                 end;
           end;
         
          B(nsol,k)=r   ; 
          w1=det(simplify(B)) ;
          u(k,1)=int(w1/w,x)          
    end;
       yp=u(1,1)*y(1);
      
       for k=2:nsol
          yp=yp+y(k)*u(k,1);
          
        end;
   disp(' The particular solution is given by the method of variation:  yp= ')
    yp = simplify(yp)
%  disp(' The homogeneous  solution is given yh=')
%    yh = simplify(yh)
%    yg=yh+yp  
      return