function y=var3(y1,y2,y3,r,x);
% -  March 7 2003(RMM)
% - Solves:   D3y+P1 D2y+ p2 Dy +p3 y = r(x)
% -  y1 y2 and y3 are known solutions
% -  A matrix form is used: A A1 A2 w=|A| W1=|A1| w2=|A2|
% -   the method of variation is used to evaluate yp.
% - usage syms x y1 y2 y3 r yy ; y1=x;y2=xLn(x);r=1;yy=var3(y1,y2,y3,r,x)
% -example y'''-6y''+11y'-6y=r(x)  has y1=exp(1x),y2=exp(2x) y3=exp(3x) find yp when r=y1
syms  w w1 w2 w3 u1 u2 u3 A A1 A2 A3 yh c1 c2 c3 yp yg
%disp(y1)
%disp(y2)
%disp(r)

A=simplify([y1 y2 y3; 
   diff(y1,x) diff(y2,x) diff(y3,x);
   diff(y1,x,2) diff(y2,x,2) diff(y3,x,2)]);
disp(A)

A1=simplify([0 y2 y3; 
   0 diff(y2,x) diff(y3,x);
   r diff(y2,x,2) diff(y3,x,2)]);
disp(A1)
A2=simplify([y1 0 y3; 
       diff(y1,x) 0 diff(y3,x);
        diff(y1,x,2) r  diff(y3,x,2)]);
disp(A2)
A3=simplify([y1 y2 0; 
    diff(y1,x) diff(y2,x) 0;
    diff(y1,x,2) diff(y2,x,2) r]);
disp(A3)

w=simplify(det(A))
w1=simplify(det(A1))
w2=simplify(det(A2))
w3=simplify(det(A3))
u1=int(w1/w,x)
u2=int(w2/w,x)
u3=int(w3/w,x)

  disp(' The particular solution is given by the method of variation')
  disp( ' yp= y1* u1(x)+y2*u2(x)+y3*u3(x) ')
  yp = simplify(y1*u1+y2*u2+y3*u3)
  
  yh=c1*y1+c2*y2+c3*y3;
  fprintf(' The Homogeneous solution yh=c1 y1+c2 y2+c3 y3 = ')
  disp(yh)
  disp(' The  general solution is given by yg=yh+yp = ')
         yg=yh+yp
   
      return