a0=input('Enter a0<>0: ');
a1=input('Enter a1<>0: ');
b0=input('Enter b0<>0: ');
discr=a1^2/4-a0;
tim=0:0.1:10;
sigma=a1/2;
if discr<0
    omega=(-discr)^0.5;
    y=b0/a0*(1-exp(-sigma*tim).*(cos(omega*tim)+sigma/omega*sin(omega*tim)));
elseif discr==0
    y=b0/a0*(1-exp(-sigma*tim).*(1+sigma*tim));
else
    root=discr^0.5;
    lambda1=-a1/2+root;
    lambda2=-a1/2-root;
    y=b0/a0*(1+lambda2/(lambda1-lambda2)*exp(lambda1*tim)+lambda1/(lambda2-lambda1)*exp(lambda2*tim));
end
plot(tim,y);
 
    