% p=0.5;
% q=6;
% omega=2;
clf
p=input('\ny"+py''+qy=cos wt \n \nEnter value of p ');
q=input('Enter value of q ');
omega=input('Enter value of omega ');
t=linspace(0,40,300);
D=p^2-4*q;
if D>0
    % two real solns
    alpha=sqrt(D)/2;
    lambda1=-p/2+alpha;
    lambda2=-p/2-alpha;
    yc=exp(lambda1*t)+exp(lambda2*t);
elseif D==0
    % repeated
    lambda=-p/2;
    yc=exp(lambda*t)+40*t.*exp(lambda*t);
else
    % complex
    beta=-p/2;
    alpha=sqrt(-D)/2;
    yc=exp(beta*t).*(cos(alpha*t)+sin(alpha*t));
end
subplot(3,1,1),plot(t,yc,'LineWidth',2)
params=[' p=',num2str(p),', q=',num2str(q)];
title(['Homogeneous solution',params],'FontSize',20)
pause

if p==0 & q==omega^2
    yp=t.*sin(omega*t)/(2*omega);
else
    A=1/sqrt((q-omega^2)^2+omega^2*p^2);
    theta=atan(-omega*p/(q-omega^2));
    yp=A*cos(omega*t+theta);
end
subplot(3,1,2),plot(t,yp,'LineWidth',2)
params=[' p=',num2str(p),', q=',num2str(q),', omega=',num2str(omega)];
title(['Particular solution',params],'FontSize',20)
pause

subplot(3,1,3), plot(t,yc+yp,'LineWidth',2)
params=[' p=',num2str(p),', q=',num2str(q),', omega=',num2str(omega)];
title(['Solution',params],'FontSize',20)