/* regular perturbations on Mathieu's equation */

 /* set up expansions */
 trunc:read("truncation order")$
 depends(x,t);
 xx:sum(x[i]*e^i,i,0,trunc);
 delta:sum(d[i]*e^i,i,0,trunc);
 /* which transition curve? */
 results:[x[0]=1,d[0]=0];
 /* prepare null function for later */
 null(t):=0;

 /* plug into d.e. and collect terms */
 de1:diff(xx,t,2)+(delta+e*cos(t))*xx$
 de2:taylor(de1,e,0,trunc)$
 for i:0 thru trunc do eq[i]:coeff(de2,e,i)$

 /* main loop */
 for i:1 thru trunc do (
 temp1:ev(eq[i],results),

 /* apply trig identities */
 temp2:expand(trigreduce(expand(temp1))),

 /* pick off secular terms */
 temp3:subst(null,cos,temp2),
 temp4:ev(temp3,x[i]=0),

 /* solve for d[i] */
 temp5:solve(temp4,d[i]),
 results:append(results,temp5),
 temp6:ev(temp2,results),
 temp7:ratsimp(temp6),

 /* use o.d.e. package to find x[i] */
 temp8:ode2(temp7,x[i],t),
 /* zap arbitrary constants */
 temp9:ev(temp8,%k1=0,%k2=0),
 results:append(results,[temp9]))$
 /* end of main loop */

 /* output results */
 deltafinal:ev(delta,results);