restart; # # Set up the exact problem # de:= (f,v)-> diff(f(v), v, v, v) + (4/v)*diff(f(v), v, v) + (2/v^2)*diff(f(v), v)-9*(4+10*v^3+3*v^6)*f(v)=0: bcs:= g(0)=1, D(g)(0)=0, (D@@2)(g)(0)=0: # # So the "exact" solution for the problem is # exact:=dsolve([de(g,x), bcs]); # # Compute what I assume to be the approximate # solution # So setup - # n:=12: lambda :=-s^3/x+s^2*(1+ln(s/x)): y := 1: # # Now compute it # for k from 0 to n do yI[k] := eval(y[k], x = s); y[k+1] := y[k] + int ( lambda*( diff(yI[k],s,s,s)+4*(diff(yI[k], s, s))/s + 2*(diff(yI[k], s))/s^2 - 9*(3*s^6+10*s^3+4)*yI[k] ), s = 0 .. x ): end do: approx := evalf( add ( y[n], i = 0 .. n ) ); # # Plot the "exact" solution in red and the "approximate" # solution in blue - can't say that the approximate # solution is particularly close - particularly for # negative x # plot( [rhs(exact), approx], x=-5..5, color=[red, blue], view=[-5..5, -20..20] ); # # Idle curiosity - what do the individual terms in # the approximate solution look like? # plot( [seq( y[j], j=0..n )], x=-5..5, view=[-5..5, -20..20] ); # # Compute the discrepancy at a couple of points - # just for interest # eval( rhs(exact)-approx, x=0.0); eval( rhs(exact)-approx, x=0.5); LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic=