function y=linode22(a,b,c,q1);
% - Revision of March 16 2002(RMM)
% - The linode22(a,b,c,q1) solves:    a D2y+ b Dy +c y =q1(x)
% - a  b c are constants and q1 is a function of x
% - the classification is used to find y1 and y2 then 
% - the method of variation is used to evaulate yp.
% - usage a=1;b=0;c=1;q1=sin(x);y=linode22(a,b,c,q1)
% -
syms x yh yp w m kar m1 m2 y1 y2 c1 c2
    disp(' Second order Linear homogeneous ODE with constant coefficients')
    disp(' a D2y + b Dy + C y= q1')  
    sprintf(' where a = %d',a)
    sprintf('       b = %d',b)      
    sprintf(' and,  c = %d',c) 
    
    disp('The characteristic equation is  a m^2+bm+c=0')
    kar=factor(a*m^2+b*m+c);
    disp('THE CHARACTERISTIC Eqn is factored AS: '); disp(kar)
      % - we  could use solve here ...    solve('a*m^2+b*m+c', 'm');
      %
      delta = b^2-4*a*c
         sprintf(' delta= b^2-4ac =%d ',delta)
   
      if delta==0  
         disp(' class B (Critical.  Double root.  m1=m2=-b/(2a)')
         m1=-b/(2*a);
         m2=m1;
         y1=exp(m1*x)
         y2=x*y1
      end
      
      if delta >0
            disp(' class A supercritical, two different real roots,')
            m1= (-b+sqrt(delta))/(2*a);
            m2= (-b-sqrt(delta))/(2*a);
            sprintf(' m1=%d',m1)
            sprintf('m2=%d',m2)
            y1=exp(m1*x)
            y2=exp(m2*x)
         end
         if delta <0
         disp(' class C subcritical , complex conjucate roots')
            alpha=-b/2/a
            beta=sqrt(-delta)/2/a
            y1=exp(alpha*x)*cos(beta*x)
            y2=exp(alpha*x)*sin(beta*x)
         end
                 yh=c1*y1+c2*y2;
         disp(' the solutions y1 and y2 are linearly independent.')
         disp(' The wronskian is given by w=y1*df(y2,x)-y2*df(y1,x)=')
         w=y1*diff(y2,x)-y2*diff(y1,x)
         w=simplify(w)
         disp(' The particular solution is given by the method of variation')
         disp( ' y= y1* u1(x)+y2*u2(x) ')
         u1=int(-y2*q1/(a*w),x)
         u2=int(y1*q1/(a*w),x)
         yp = simplify(y1*u1+y2*u2)
         fprintf(' The particular solution')
         
         fprintf(' The Homogeneous solution yh=c1 y1+c2 y2 =')
         disp(yh)
         
         disp('The  general solution is given by y=yh+yp =')
         y=yh+yp;
      return