I have to integrate numerically the next system of partial differential equations, in "r" and "z" coordinates (steady state, cylindrical coordinates): -Continuity equation -Momentum equation (with viscous and convection terms) -Energy equation (with heat transfer, convection, and internal generation terms). I heard something about Crank-Nicholson method. I have to implement this problem in Matlab or Fortran. Anybody know about how can I solve this heavy task? Anybody know about some web where I could obtain a method for solving it, and implemented in the software specified?. Thanks for answering.