function cc4plot(a,b,c,d,e,y0,yp0,ypp0,yppp0,t0,tend,g)
% cc4plot(a,b,c,d,e,y0,yp0,ypp0,yppp0,t0,tend,g) -- Plot the solution to
% the fourth order linear constant coefficient ODE 
% ay^(4) +by''' + cy'' +dy' +ey =g(t),
% with 
% y(t0)=y0, y'(t0)=yp0, y''(t0)=ypp0, y'''(t0)=yppp0 
% from t=t0 upto t=tend. 
% (The function g is optional, and is '0' by default)
%
% Examples:
%         Fig4.2.2: cc4plot(1,0,0,0,-1,7/2,-4,5/2,-2,0,15)
%         Fig4.2.3: cc4plot(1,0,0,0,-1,7/2,-4,5/2,-1.99999,0,15)
  
  if nargin < 12
    g='0';
  end
  
  set(0,'DefaultAxesLineWidth',3);
  set(0,'DefaultLineLineWidth',2);
  get(0,'Default');
  options = odeset('OutputFcn',@odeplot,'OutputSel',1);
  ode45(@(t,y) ode(t,y,a,b,c,d,e,g) ,[t0 tend],[y0 yp0 ypp0 yppp0],options);
  set(gca,'FontSize',18)
  xlabel('t','FontSize',18)
  ylabel('y','FontSize',18)
  r=roots([a b c d e]);
  r1=num2str(r(1)); r2=num2str(r(2)); r3=num2str(r(3)); r4=num2str(r(4));
  title(sprintf(['%g y"" + %g y"'' + %g y" + %g y'' + %g y = %s\nroots=%s,%s,%s,%s,',...
		 '\ny(0)=%g,y''(0)=%g,y''''(0)=%g,y''''''(0)=%g'],...
		a,b,c,d,e,g,r1,r2,r3,r4,y0,yp0,ypp0,yppp0),'FontSize',18);
  
function dydt=ode(t,y,a,b,c,d,e,g)

  gf=inline(g);
  gv=feval(gf,t);
  dydt = [y(2);y(3);y(4);(gv-b*y(4)-c*y(3)-d*y(2)-e*y(1))/a];
