How to speed-up a code containing several symbolic integrations and derivatives?

edit

Hi,

After days of confusion now I have a sage code ready to be used. It consists of a set of iterative equations in a loop, at each step the equations --containing a number of symbolic integration and differentiation-- act on a previous estimation of the unknowns to update them to a newer estimation.

Sounds good up to this point. However, depending on the initial estimation that I introduce to the iterative loop the code may run fast (like when I give a constant or simply x as initial estimation) or too slow (like when I introduce sin(x) or other even simple polynomials). I understand that having such integrations and differentiations make the expressions rapidly grow large so that the computations would be time consuming but it is extremely slow, not one minute or 10 minutes or so, I just waited for some hours and sage did not answered anything new, and I doubt if the 100% cpu usage shown by the system's monitoring app is really spent on the calculations and sage can give the answer at the end.

As much as I read from documentations and elsewhere adding a %cython command at the beginning of the code can speed up the code but it gives me back an error e.g. of the form "undeclared name not builtin: sin". Does %cython command work also with symbolic calculations? Also I saw something about "fast_callable" function and sympyx, what about they, if they can help speed up such calculations where can I find a good intro toward them?

Best regards

EDIT.

Here is the code:

var('x,y,z,t, x1,y1,z1,t1, x2,y2,z2,t2, x3,y3,z3,t3') 

N=5         # N: number of iterations
            # or N=int(input("How many iteration do you need for your Adomian Decomposition code to run? "))
Re=1000000  # Re: global Reynolds number
x_B1=-1; x_B2=1; y_B1=-1; y_B2=1; z_B1=-1; z_B2=1       # the non-dimensional boundaries of the computational region!
T=10        # T: the ending non-dimensional time of computations
# p0, v0: used in definition of initial condition for the unknown functions

assume(x_B1<=x, 0<t)     # needed for our integrations from X_B1 to x and from 0 to t !


#------------------------------------------------
# the following lists are defined for, 1st, easier addressing and, 2nd, more compact importing
# of terms like "sum_j u_j du_i/dx_j" as otherwise x_j was not known to the code
var('q')                # q is a very dummy variable to fill the zeroth place in all my lists to avoid confusion while coding, instead of always using i-1 or j+1 as indices!
R0=[q,x,y,z,t]          # Only R0 needs a zeroth element of q
R1=[x1,y1,z1,t1]
R2=[x2,y2,z2,t2]
R3=[x3,y3,z3,t3]


#--------------------------
# empty lists that will be filled with functions defined in the loop below
phi0=[q]
phi1=[q]
phi2=[q]
#-----------------
for i in range(1,5):        # i = 1,2,3,4  / indices 1,2,3 are reserved for velocity ...