@vv 

Thank you, 
In the system of equations A , B and C must be replaced by A(t), B(t) and C(t) 
because they unknown functions 

I run the same code I get 
Error, (in ODEtools/ODESolStruc/info) incorrect ODESolStruc: expected 1 old dependent variable in the 2nd set of transformation equations; received: 0
code_solve_system_of_equations_vv.mw

I have 8 unknows functions, \Zeta_i(t,x,y,z), i=0,1,2,3 and A(t), B(t), C(t), \psi(t,x,y,z)

Thank you for your help 

@C_R 

The difference equation is the first equation in the maple code : system defined at three different times (n+1), n and n-1 

@vv 

thank you

@C_R 

I think everything is well now, in my code. 

A_less_than_one.mw

In the definition of beta_1 the loop from 1 to q is correct, beta_1 is a coefficient used  in A 

and everything run correctly 

thanks

@Carl Love 

so only one i a field, so can not get a homomorphism between them

@tomleslie 

This means that in geneal case, for a given x,y,z and w , we can always find n, a,b and c solution of the system 

Since we have many choice of x , y, z and we so for each choice we have a solution 

there are infintely many solution of the system 

@Axel Vogt 

But Maple return a zero solution 
Instead of a nonzero solution 

for a fixed x=1, y=1, z=6 and w=-1 we have 

n=4, a=-1, b-0, c=5 

is a solution for this fixed vector defined by x,y,z and w 

@vv 

Big thanks 

@vv 

The function f is defined only in the unit ball centered at zero and outside is zero. 
So, we have singularity at zero, and it is integrable in the vicinity of zero since 2< 3 

@vv 

Using polar coordinate you will integrate over the range [0,1]  the function r^(-2) r^(3-1) =1  and so integrable 

In general case, f(x)=|x|^(-alpha) times indicator of unit ball, is integrable over the unit ball if alpha< dimension of space  ( in our example alpha equal 2 < 3 ), so the function is integrable 

@Zeineb 

Psi(x,y), theta(x,y) converge to constant constant when y goes to infinity 

@Zeineb 

All my function are continuous, and the remaining are positve parameter 

@vv 

Thank you for reading my code, 

Pe is not defined as a function, Pe is a real number, we can replace Pe by Z , 

@dharr 

Thnak you for the idea. 
Using Routh-Hurwitz 'array'  can we solve the system to get the set of values where the roots are in the left half plane.

@mmcdara 

Thanks for your remarks. 

All parameters are positive in my equation. But, Maybe we can add  a condition so that we can determine or know the sign of the roots of the equation. 

My first goal is to study the stability of a given equilibrium point, so I computed the characteristic equation which contains  many paramters. So, if I know the sign of my roots I can know the nature of my equilibrium point. 

Maybe we can add some condition that gives the sign of the roots