from sympy import *
U,V,x,y = symbols('U,V,x,y')
L = symbols('L')
init_printing()
get_ipython().magic(u'matplotlib inline')
# Blasius equation
beq = 2*Derivative(U(y),y,3)+Derivative(U(y),y,2)*U(y)
Eq(beq,0)
# ## substitution of variables
# from y to x
#new variable: x
#old variable: y
y2x = 2*y/(L+y)-1
# from x to y
x2y = solve(y2x-x,y)[0]
# ## Dictionary for change of variable
sd = zip([diff(U(y),y,i) for i in range(4)],[diff(y+V(y2x),y,i).subs({y:x2y}).simplify().doit() for i in range(4)])
sd
# Change of variables from high order derivative to low order ones.
#
# *sympy* function **replace** is used rather than **subs** which will give incorrect result.
for sdi in sd[::-1]:
beq = beq.replace(*sdi)
beq=beq.simplify()
# ## some hacks to generate code for chebfun.chebop in matlab
DIFF=symbols('diff')
sd2 = zip([Derivative(V(x),x,i)for i in range(4)],[DIFF(V,i) for i in range(4)])
beq2 = beq.args[-1]
for sdi in sd2[::-1]:
beq2 = beq2.replace(*sdi)
beq2
print octave_code(beq2.subs({L:1}))