/*

(d12) This program computes a perturbation solution for the


limit cycle in Van der Pol's equation.  Call it by typing:


                     VANDERPOL(N)


where N is the order of truncation.
*/


vanderpol(n):=(setup1(n),setup2(n),
	  for i thru n do
	      block(step1(i),step2(i),IF i > 1 THEN output1(i),
		    IF i = n THEN go(end),step3(i),step4(i),step5(i),end),
	  output2(n))$
setup1(n):=(w:1,for i thru n do w:w+k[i]*e^i,x:2*cos(t),
       for i thru n do x:x+y[i](t)*e^i)$
setup2(n):=(temp1:diff(x,t,2)+x/w^2-e*(1-x^2)*diff(x,t)/w,
       temp1:taylor(temp1,e,0,n),for i thru n do eq[i]:coeff(temp1,e,i))$
step1(i):=temp1
      :expand(trigreduce(expand(ev(eq[i],makelist([e[j],f[j]],j,1,i-1),
				   diff))))$
step2(i):=(
      IF i = 1 THEN f[i]:solve(coeff(temp1,cos(t)),k[i])
	  ELSE f[i]:solve([coeff(temp1,cos(t)),coeff(temp1,sin(t))],
			  [k[i],b[i-1]]),temp1:ev(temp1,f[i]))$
step3(i):=temp1:ev(ode2(temp1,y[i](t),t),%k1:a[i],%k2:b[i])$
step4(i):=(temp2:rhs(temp1),temp2:diff(temp2,t),
      temp2:solve(ev(temp2,t:0),a[i]))$
step5(i):=e[i]:ev(temp1,temp2)$
output1(i):=(print(expand(ev(e[i-1],f[i]))),print(" "))$
output2(n):=(print("w=",ev(w,makelist([f[j]],j,1,n))),print(" "))$