function y=var4(y1,y2,y3,y4,r,x);
% -  March 9 2003(RMM) Fourth order linear ode  by variation.
% -   Solves:   D4y+P1 D3y+ p2 D2y +p1Dy+ Qy= r(x)
% -  y1, y2, y3 and y4 are known solutions linearly independent solutions
% -  A matrix form is used: A A1 A2 A3 w=|A| W1=|A1| w2=|A2| w3=|A3| w4=|A4|
% -   the method of variation is used to evaluate yp.
% - usage syms x y1 y2 y3 y4 r y ; y1=x;y2=xLn(x);y3=1;y4=exp(x);r=1;y=var4(y1,y2,y3,y4,r,x)
% -example y''''-6y'''+11y''-6y'=r(x)  has y1=exp(1x),y2=exp(2x) y3=exp(3x) y4=1 find yp when r=y1
syms  w w1 w2 w3 w4 u1 u2 u3 u4 A A1 A2 A3 A4 c1 c2 c3 c4 yh yp yg
%disp(y1)
%disp(y2)
%disp(y3)
%disp(y4)
%disp(r)

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

A1=simplify([        0             y2           y3           y4; 
                     0        diff(y2,x)   diff(y3,x)   diff(y4,x);
                     0        diff(y2,x,2) diff(y3,x,2) diff(y4,x,2);
                     r        diff(y2,x,3) diff(y3,x,3) diff(y4,x,3)]);
disp(A1)

A2=simplify([        y1        0          y3           y4; 
                diff(y1,x)     0     diff(y3,x)   diff(y4,x);
                diff(y1,x,2)   0     diff(y3,x,2) diff(y4,x,2);
                diff(y1,x,3)   r     diff(y3,x,3) diff(y4,x,3)]);
disp(A2)

A3=simplify([        y1         y2              0         y4; 
                diff(y1,x)   diff(y2,x)         0    diff(y4,x);
                diff(y1,x,2) diff(y2,x,2)       0    diff(y4,x,2);
                diff(y1,x,3) diff(y2,x,3)       r    diff(y4,x,3)]);
disp(A3)
A4=simplify([       y1         y2              y3           0; 
                diff(y1,x)   diff(y2,x)   diff(y3,x)        0;
                diff(y1,x,2) diff(y2,x,2) diff(y3,x,2)      0;
                diff(y1,x,3) diff(y2,x,3) diff(y3,x,3)      r]);
disp(A4)
w=simplify(det(A))
w1=simplify(det(A1))
w2=simplify(det(A2))
w3=simplify(det(A3))
w4=simplify(det(A4))
u1=int(w1/w,x)
u2=int(w2/w,x)
u3=int(w3/w,x)
u4=int(w4/w,x)
  disp(' The particular solution is given by the method of variation')
  disp( ' yp= y1* u1(x)+y2*u2(x)+y3*u3(x)+y4*u4(x) ')
  yp = simplify(y1*u1+y2*u2+y3*u3+y4*u4)
  yh=c1*y1+c2*y2+c3*y3+c4*y4;
  fprintf(' The Homogeneous solution yh=c1 y1+c2 y2+c3 y3+ c4 y4 = ')
  disp(yh)
  disp(' The  general solution is given by yg=yh+yp = ')
         yg=yh+yp
      return